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In  this  dissertation  we  consider  special  cases  of  optimal 
differential  problems,  the  dual  of  the  uncapacitated  minimum  cost 
network  flow  problem  and  project  scheduling  problem,  and  develop 
algorithms  for  solving  these  problems. 

We  develop  three  dual  simplex  algorithms  for  the  minimum  cost 
flow  problem,  and  its  special  cases,  the  transportation  problem  and  the 
assignment  problem.  The  first  algorithm  provides  a polynomial  bound  for 
the  assignment  problem.  The  other  two  algorithms  use  scaling  techniques 
and  provide  a polynomial  bound  for  the  minimum  cost  flow  problem  and 
the  transportation  problem. 
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We  study  the  project  scheduling  problem  with  the  objective  of 
maximizing  the  present  value  of  the  cash  flows  and  develop  an  algorithm 
which  solves  this  nonlinear  programming  problem  optimally. 

We  present  computational  experience  obtained  by  comparing  the 
algorithms  with  related  algorithms  in  the  literature. 
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CHAPTER  1 
INTRODUCTION 


1.1  Purpose  of  the  Study 

There  are  a number  of  special  structured  mathematical  programming 
models  that  are  referred  to  as  "network  models."  These  special 
structures  can  be  exploited  to  allow  the  construction  of  efficient 
algorithms.  In  this  research,  we  will  examine  some  network  models  that 
are  referred  to  as  optimal  differential  problems  [1] . 

Optimal  differential  problems  are  those  optimization  problems  in 
which  the  contraint  set  consists  of  differentials  of  the  potentials. 
Rockafellar  [1]  discusses  some  special  cases  of  linear  differential 
problems.  These  special  cases  include  the  duals  of  transportation 
problems,  the  duals  of  assignment  problems,  and  project  scheduling 
problems.  He  presents  a simplex  algorithm  for  solving  these  problems. 
For  the  special  cases,  this  algorithm  is  equivalent  to  the  dual  simplex 
algorithm.  However,  a specific  rule  for  pivoting  is  not  given  for  this 
simplex  algorithm. 

In  this  study,  we  will  consider  the  dual  simplex  approach  in 
solving  the  uncapacitated  minimum  cost  flow  problem  with  its  special 
cases,  the  transportation  problem  and  the  assignment  problem.  Recently, 
dual  simplex  approaches  in  solving  these  problems  have  been  studied  in 
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the  literature  in  the  search  for  finding  polynomially  bounded 
algorithms . A particular  leaving  arc  selection  rule  allowed  us  to 
develop  dual  simplex  algorithms  with  polynomial  bounds.  Moreover,  an 
algorithm  that  solves  the  project  scheduling  problem  with  the  objective 
function  of  maximizing  the  present  value  of  the  cash  flows  is 
developed.  This  algorithm  uses  ideas  similar  to  that  of  the  dual 
simplex  algorithm  for  the  minimum  cost  flow  problem. 

1.2  Scone  of  the  Study 

In  each  of  the  remaining  chapters,  we  will  concentrate  on  one 
specific  problem.  Each  problem  will  be  defined,  and  the  related 
literature  will  be  discussed  within  the  chapter.  Computational 
experimentation  comparing  the  discussed  algorithms  with  related 
algorithms  in  the  literature  will  be  reported. 

In  Chapter  2,  we  will  develop  three  polynomially  bounded  dual 
simplex  algorithms  for  the  minimum  cost  flow  problem.  In  two  of 
algorithms  the  scaling  approach  will  be  considered. 

In  Chapters  3 and  4,  special  cases  of  the  minimum  cost  flow 
problem,  transportation  and  assignment  problems,  will  be  studied. 

In  Chapter  5,  the  project  scheduling  problem  is  studied,  and  an 
algorithm  which  solves  this  nonlinear  programming  problem  optimally  is 
presented. 

The  notation  used  throughout  the  dissertation  will  be  summarized 


in  the  Appendix. 


CHAPTER  2 

THE  MINIMUM  COST  FLOW  PROBLEM 


The  minimum  cost  flow  problem  arises  whenever  items  must  be 
shipped  from  several  sources  through  a network  in  response  to  demand 
with  the  minimum  cost  of  shipment.  Problems  of  this  type  occur  in  the 
analyses  of  several  systems,  including  production-distribution  systems, 
urban  traffic  systems,  railway  systems,  military  logistic  systems  and 
communication  systems. 

2,1  Formulation  of  the  Problem 

Consider  a network  representation  g' ,e' ) of  a minimum  cost 
flow  problem  with  V =Vg  u Vp,  where  Vg  is  the  set  of  supply  nodes  and 

I 9 / 

Vp  is  the  set  of  demand  nodes.  If  there  exists  a direct  arc  from  ieV 
to  jev  , then  (i,j)eE  . For  each  arc  (i,j)eE  , we  associate  an  integer 
cjj  which  corresponds  to  a unit  shipping  cost  over  the  arc  (i,j). 
Moreover,  with  each  node  i€Vg,  we  associate  a positive  integer  a ^ 
designating  the  supply  at  that  point  and  with  each  node  jeVp,  we 
associate  a nonnegative  integer  bj  corresponding  to  the  demand  at  point 
j.  If  bj  =0,  then  node  j is  called  a transshipment  node. 

For  convenience,  for  each  node  ieV*  , let  oc(i)  represent  the 
"forward  star"  of  node  i which  is  the  set  of  all  the  outgoing  arcs  from 
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node  i,  and  let  /3(i)  represent  the  "backward  star"  of  node  i which  is 
the  set  of  all  the  incoming  arcs  to  node  i. 


The  minimum  cost  flow  problem  (MCP)  can  now  be  stated  as  a linear 
programming  problem  : 

(MCP)  Minimize  2 (2.1) 

(i.j)eE' 

subject  to 


a^  if  ieVg 


2 x^j  - 2 Xjj 

(i, j)ex(i)  (i, j)e£(i) 


(2.2) 


-b-^  if  ieVj) 


xij  * 0 


(i,j)eE  . 


(2.3) 


In  what  follows  we  assume  that  the  sum  of  supplies  is  equal  to 
the  sum  of  demands,  i.e., 

2 ^ai  - 2 (bj  . (2.4) 

ieVs  jeVD 

In  Section  2.3,  we  present  the  dual  of  (MCP).  In  Section  2.4,  we 
provide  some  definitions  on  the  tree  structure  of  a dual  feasible 
solution.  Also  in  this  section,  we  present  the  idea  of  a dual  simplex 
pivot,  the  "induced  subtree"  of  a given  node  and  a general  overview  of 
the  algorithm  and  present  the  algorithm  for  solving  the  minimum  cost 
flow  problem.  Section  2.4  also  presents  the  validity  and  the 
computational  bound  of  the  algorithm.  In  Sections  2.5  and  2.6,  two 
scaling  algorithms  with  polynomial  bounds  are  developed  for  the  minimum 
cost  flow  problem. 
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2.2  Literature  Review 

The  transportation  problem,  a special  case  of  minimum  cost  flow 
problem,  was  first  solved  by  Hitchcock  [2]  in  1941  and  Koopmans  [3]  in 
1946.  In  1951,  Dantzig  [4]  developed  a variant  of  the  simplex  algorithm 
for  the  transportation  problem.  Orden  [5]  extended  these  results  for 
the  minimum  cost  flow  problem. 

Early  approaches  other  than  the  primal  simplex  that  have  been 
proposed  are  out-of -kilter  by  Fulkerson  [6],  primal-dual  by  Ford  and 
Fulkerson  [7],  dual  by  Balas  and  Hammer  [8],  flow  augmentation  by 
Busacker  and  Gowen  [9],  and  negative  cycle  method  by  Klein  [10], 

In  1950s  and  early  1960s,  the  out-of -kilter  algorithm  was  seen  as 
superior  to  all  the  other  methods.  Development  of  fast  and  efficient 
codes  for  the  primal  problem  begun  in  1970 's  with  the  development  of 
new  computer  science  tools.  Currently,  based  on  the  experiments  which 
compare  the  different  algorithms  for  minimum  cost  flow  problems , it  is 
believed  that  primal  simplex  algorithms  are  superior  to  all  the  others 
[11,12]  since  primal  implementations  are  faster  and  require  less 
storage . 

The  successful  implementations  of  the  primal  simplex  network 
algorithms  are  based  on  compact  representation  of  the  basis,  efficient 
ways  of  determining  the  entering  and  leaving  arcs  and  efficient 
techniques  to  update  the  simplex  multipliers  at  each  pivot. 

Different  data  structures  have  been  used  to  represent  the  basis. 
In  addition  to  predecessor  and  thread  labels,  the  following  labels  are 
proposed;  "number  of  nodes"  label  by  Glover  and  Klingman  [13], 
"augmented  predecessor"  label  by  Glover  , Karney  and  Klingman  [14]  , 
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"augmented  thread"  label  by  Glover,  Klingman  and  Stutz  [15],  "preorder 
distance"  label  by  Bradley,  Brown  and  Graves  [16].  Ali  et  al.[17] 
conclude  that  data  structures  which  use  the  thread  and  predecessor 
labels  in  conjunction  with  other  labels  are  the  most  efficient  data 
structures  available  for  simplex  specializations. 

Efficient  methods  to  perform  pricing  for  determining  the  entering 
arc  generally  include  a candidate  queue  [16,18]  which  is  used  to  price 
only  a subset  of  the  nonbasic  variables  at  each  pivot. 

Even  though  primal  simplex  algorithms  are  superior  to  other 
methods  computationally,  they  may  fail  to  terminate  due  to  an  infinite 
number  of  degenerate  pivots  caused  by  cycling.  In  order  to  prevent 
cycling,  Cunningham  [19] , and  Barr,  Glover  and  Klingman  [20]  proposed 
algorithms  which  work  only  with  strongly  feasible  trees.  A strongly 
feasible  tree  is  defined  to  be  a feasible  tree  in  which  every  arc  with 
a flow  of  zero  is  directed  away  from  the  root.  Cycling  is  prevented  by 
a specific  choice  of  the  leaving  arc.  Bland's  [21]  rule  for  choosing 
the  entering  and  leaving  arcs  (choosing  the  variable  with  the  least 
index  possible)  also  prevents  cycling. 

In  the  absence  of  cycling,  primal  simplex  algorithms  can  still 
encounter  an  exponentially  long  sequence  of  consecutive  degenerate 
pivots.  This  phenomenon  is  called  stalling.  In  order  to  avoid  stalling, 
Cunningham  [22]  proposed  to  use  some  entering-arc  rules  combined  with 
the  maintenance  of  strongly  feasible  trees.  For  this,  each  arc  is 
considered  periodically  as  a choice  for  entering  arc.  The  rules  that 
maintain  these  are  Least  Recently  Considered,  Least  Recently  Basic,  and 
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Inward  Most  Negative  rules.  Rules  similar  to  these  are  also  suggested 
in  other  studies  [12,23,24]. 

However,  Zadeh  [25]  proved  that  when  Dantzig's  rule  for  the 
entering  variable  is  used,  the  primal  simplex  algorithm  may  take  an 
exponential  number  of  pivots  on  some  problems  while  all  the  bases 
remain  nondegenerate.  Also  he  presented  bad  networks  for  primal-dual, 
path  and  cycle  methods. 

The  algorithms  that  are  of  primal-dual  type  include  the 
following : 

Tomizawa  [26]  developed  a variant  of  the  ordinary  primal-dual 
algorithm  in  which  the  associated  shortest  path  problem  was  solved  by 
Dijkstra's  algorithm. 

Bertsekas  [27]  introduced  a unified  framework  for  primal-dual 
methods  in  minimum  cost  network  flow  problems  by  giving  different 
procedures  for  implementing  flow  augmentations,  price  adjustments  and 
ascent  steps. 

Following  Zadeh' s findings,  algorithms  other  than  primal  simplex 
and  primal-dual  have  been  considered  in  the  hopes  of  finding 
polynomially  bounded  algorithms . 

Edmonds  and  Karp  [28]  were  the  first  to  solve  minimum  cost 
network  flow  problems  in  polynomial  time.  Their  algorithm  uses  a 
scaling  technique  and  solves  a sequence  of  0[r]  different  network  flow 
problems  where  2r  < c(u,v)  for  all  arcs  (u,v)  and  c(.,.)  is  the 
capacity  of  an  arc.  These  problems  differ  in  their  arc  capacities  where 
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for  problem  p,  capacities  are  obtained  by  deleting  the  p low-order 
digits  in  the  binary  representations  of  the  original  capacities.  The 
algorithm  has  a computational  bound  of  0(r|V|^).  Orlin  [29]  proved  that 
a variant  of  the  Edmonds -Karp  scaling  algorithm  solves  the  minimum  cost 
network  flow  problem  in  0( ( | V | ^log | V | ) ( | E | + | V | log | V | ) , a bound 
independent  of  input  data. 

Rock  [30]  described  a variant  of  the  Edmonds  and  Karp  scaling 
algorithm  with  a computational  bound  of  0(oc|V|^)  where  2“  < C(u,v)  for 
all  arcs  (u,v)  and  C(.,.)  is  the  cost  of  an  arc.  Here  the  bound  depends 
on  the  size  of  the  costs. 

Orlin  [29]  also  developed  two  genuinely  polynomial  dual  simplex 
algorithms  for  the  minimum  cost  flow  problems  that  have  bounds  of  0[U 
(|V|log  | V | + | E | ) ] and  0[U*(|V|)]  pivots  where 

U*  = min  [|V|2  log  |V| , |V|BIT(b)].  A genuinely  polynomial  algorithm  is 

an  algorithm  that  has  a computational  bound  which  is  polynomial  in  the 
number  of  nodes  and  is  independent  of  both  costs  and  capacities.  Here 
BIT(b)  is  the  minimum  integral  value  of  s such  that  2sb  is 
integer-valued.  In  Orlin' s formulation  of  the  problem,  the  vector  b is 
the  vector  of  weights  of  the  nodes  and  it  is  scaled  so  that 
| b £ | < 1.  The  computational  bounds  of  these  algorithms  are  0(U*|V|^) 
and  0(U  |Vp)  respectively.  Orlin' s algorithms  use  the  idea  of  the 
Edmonds  and  Karp  scaling  method. 

Armstrong,  Klingman  and  Whitman  [31]  introduced  a variant  of  the 
dual  simplex  algorithm  in  which  different  list  structures  for  storing 
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the  basis  were  considered  and  multiple  pivot  strategy  was  employed.  No 
computational  bounds  were  given  for  their  algorithm. 

Ikura  and  Nemhauser  [32,33]  provided  a dual  simplex  algorithm  for 
the  transportation  problem  with  an  upper  bound  on  the  simplex  pivots 
polynomial  in  the  length  of  the  input  data.  More  precisely,  the 
algorithm  has  a bound  of  0(|V|2M)  simplex  pivots  where  M is  the  stun  of 
the  supply  and  demand.  Moreover,  if  the  Edmonds  and  Karp  algorithm  is 
used  in  a recursive  manner  to  adjust  the  supplies  and  demands,  the 
upper  bound  on  the  number  of  simplex  pivots  of  the  revised  algorithm 
becomes  a polynomial  of  0(|V|^(r+l))  where  r is  defined  as 
2T  < max  wj  < 2r+1.  Here  wj  is  the  weight  (supply  or  demand  ) of 
node  j.  This  algorithm  has  a computational  bound  of  0( | V | 5 (r+1) ) . 

Recently,  Tardos  [34]  develeloped  a genuinely  polynomial 
algorithm  for  the  minimum-cost  circulation  problem.  Her  algorithm  has  a 
computational  bound  of  0( j E | 2T( | E | , | V j ) log  |E|).  Here  T(|E|,|V|)  is  the 
order  of  solving  a maximum  flow  problem  in  a network  with  jEj  arcs  and 
| v | nodes.  Her  algorithm  solves  at  most  jEj  subproblems  obtained  by 
scaling  the  cost  coefficients. 

Fujishige  [35]  provided  an  algorithm  for  the  dual  of  Tardos 's 
problem  which  runs  in  0(  |E]  ( |V|+m0)2log  jEj)  time  where  m0  is  the 
number  of  arcs  with  both  finite  lower  and  upper  capacities.  His 
algorithm  solves  at  most  j E | subproblems  obtained  by  scaling  the 
capacities . 
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2,3.  The  Dual  Formulation 


Consider  the  augmented  graph  G=(V,E)  which  is  obtained  from 

9 9 9 

G - (V  , E ) by  appending  a "super  source”  with  node  number  0 connected 
to  each  source  node  jeVg  via  an  arc  (0,j)  with  associated  costs 
cq j = 0.  Let  ag  - 0.  This  situation  is  depicted  in  Figure  2.1  below. 
The  types  of  nodes,  whether  a source  node  or  a demand  node,  are 
indicated  by  the  letters  S or  D before  the  node  numbers. 


G-(V.E) 

Figure  2.1:  The  graphs  G and  G 


Consider  the  dual  of  the  LP  formulation  of  (MCP)  by  using  G. 
(DMCP)  Max  E a^w,-  - E b-jW-j  (2.5) 


Max  E aiwi 

ievs 


S bjWj 

jeVD 


subject  to 


Wj  < Cij 


(i.j)eE. 


(2.6) 


11 


2.4,  A Dual  Simplex  Algorithm  for  the  Minimum  Cost  Flow  Problem 

In  this  section,  we  will  first  discuss  separately  the  main  steps 
of  the  algorithm;  initialization,  determination  of  the  leaving  and 
entering  arcs.  Then,  we  will  present  the  algorithm  in  Section  2.4.6. 

2.4.1  Notation 

Throughout  Section  2.4  we  will  use  the  following  notation. 

A path  in  G is  an  alternating  sequence  Pk=v0 > el , • • . , e^, v^  of 
nodes  and  arcs  such  that  e^  = (v^.p.vp)  or  e^  = (v^,vp_p).  In  the 
former  case,  e^  is  called  a forward  arc  of  P.  In  the  latter  case  it 

is  called  a backward  arc  of  P. 

Let  T-(V,Et)  be  a spanning  tree  rooted  at  node  0.  Define  PR(x)  to 
be  the  predecessor  function  indicating  the  node  which  precedes  node  x 
on  the  unique  path  from  node  x€V  to  node  0 in  T.  Let  D(i)  represent  a 
depth  function  defined  as  the  number  of  nodes  traversed  from  node  0 to 
node  i on  a unique  path  in  T. 

Let  Vp  be  the  set  of  nodes  ieV  such  that  the  arc  (PR(i),i)  is  a 
forward  arc  in  T.  Similarly,  let  Vg  be  the  set  of  nodes  ieV  such  that 
the  arc  (PR(i),i)  is  a backward  arc  in  T. 

2.4.2  Construction  of  the  Initial  Tree 

Let  T0=  (V,Ep)  be  a shortest  path  tree  of  G rooted  at  node  0, 
using  costs  as  the  weights.  If  all  the  costs  are  nonnegative,  we  can 
solve  the  shortest  path  problem  using  Dijkstra's  algorithm  in  0(|V|^) 
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time.  If  the  costs  are  not  nonnegative,  then  we  can  solve  the  shortest 
path  problem  using  the  label-correcting  algorithm  for  shortest  paths  in 
0(|V|3)  time.  In  this  study,  we  assumed  that  there  exists  no  negative 
cycles  in  G.  Once  we  solved  the  shortest  path  problem,  we  let  w^ 
represent  the  negative  of  the  shortest  path  from  node  0 to  node  i in 

TO- 

Proposition  1:  W^,  ieV,  represents  a basic  feasible  solution  for 
(DMCP) , the  dual  problem. 

Lemma  2.1:  If  the  tree  Tq  possesses  a feasible  solution  for  the 

(MCP)  by  only  allowing  the  tree  arcs  to  permit  flow,  then  the  solution 

is  optimal  for  the  (MCP) . 

Proof : Primal  feasibility,  dual  feasibility  and  complementary 
slackness  are  all  satisfied. 

2,4,3  Checking  Primal  Feasibility 

Let  T=(V,E'p)  be  a spanning  tree  of  G which  possesses  a dual 
feasible  solution.  We  call  a tree  T^-CV^.E^)  the  subtree  "induced"  by 
node  ieVp  if  this  subtree  is  obtained  by  dropping  arc  (PR(i),i)  and 
taking  the  subtree  that  does  not  contain  the  root  node.  (See  Figure 
2.2) 

Let  = V^g  u V^j)  where  V^g  is  the  set  of  supply  nodes  in  and 

V^j)  is  the  set  of  demand  nodes  in  T^.  Let 

™ 2 aj  - 2 bj  . 

JeViS  JeViD 


Si 


(2.7) 


13 


Lemma  2.2:  Given  a rooted  spanning  tree  T-(V,ET)  which  possesses 
a dual  feasible  solution,  if  there  exists  a node  ieVp  such  that  S^X), 
then  the  corresponding  primal  solution  is  infeasible. 

Proof . Since  S^>0,  if  arc  (PR(i),i)  is  the  dropped  arc  in 
obtaining  the  induced  subtree  T^,  then  in  the  computation  of  flows,  by 
only  using  the  tree  arcs,  we  will  have  X(pR^)  ^ - -Sp  < 0. 

2,4,4  Determination  of  the  Leaving  Arc 

Given  an  induced  subtree  Tp  — (Vp,Ep),  we  define  a subtree 
Ti~  (Vp,Ep)  of  Tp  and  its  corresponding  Sp  value  by  the  procedure 
PRUNE1.  Before  stating  the  exact  steps,  we  will  give  the  general  idea 
of  the  procedure.  Initially,  given  an  induced  subtree  Tp , we  set  Sj  to 
be  equal  to  Sj  for  all  nodes  jeTpDVp.  Then  we  start  from  the  tip  nodes 
in  Tp,  and  eliminate  all  subtrees  Tj  with  positive  Sj  values  from 
Once  we  eliminate  a subtree,  we  update  the  value  in  the  remaining 
subtree.  When  all  the  subtrees  with  positive  updated  sj  values  are 
eliminated  from  Tp,  the  resulting  tree  is  called  Tp  with  its 
associated  value. 

PROCEDURE  " PRUNE 1" 

INPUT  : Ti  - (Vi.Ei),  D(k) , keVp  n Vp. 

OUTPUT  : A subtree  T-=  (v[,e[)  c Ti  with  the  property  that  Sj  < 0 for 
each  jeV^pnVp  and  its  Sp  value. 


14 


BEGIN 

INITIALIZE  : Set  Tj  < Tj , and  Sj-Sj  for  all  jeV^nVp. 

Let  R - { j + i e n VF  : Sj  > 0 } . 

Z^=0 . 

WHILE  R + <t>  DO 

Find  jeR  such  that  D(j)  > D(k)  , for  all  kfjeR.  Obtain 

Tr  ■ 

Set  T{  - (V-\Vj,  e[\Ej  U ( (PR(j ) , j ) } )• 

Zi  <■ • • Zi  U { j } . 

Update  Sj , jeV^  n Vp  and  R. 

ENDWHILE 

Compute  S^-Saj-Ebj  (2.8) 

JeViS  JeViD 

where  and  V^q  are  the  sets  of  supply  and  demand 

nodes  in  respectively. 

END 


Note  that  Z^  = {j  : jeVp  n Vj_,  j=fi,  S j >0 } is  the  set  of  nodes 
whose  subtrees  are  pruned  from  subtree  with  the  procedure  PRUNE1. 

We  will  also  define  a set  where 


Zi  U { i] 

Zi 


if  Sj>0 


otherwise . 


Note  that  for  a node  ieVp,  Si  is  the  negative  of  the  flow  on  arc 
(PR(i),i)  which  defines  the  induced  subtree  Ti;  and  Si  is  the  negative 
of  the  flow  on  arc  (PR(i),i)  when  Ti  is  used  for  computation  of  flows. 
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This  might  be  viewed  as  the  individual  contribution  of  the  subtree 
to  the  flow  on  arc  (PR(i),i). 

Example  2,1.  Consider  the  minimum  cost  flow  problem  in  Figure 
2.2.  It  is  given  that  a^  - 10,  a2  = 20,  a3  - 20,  a4  = 10,  and  b5  = 5, 
t>6  = 3,  bj  = 1,  bg  - 1,  bg  - 0,  b3Q  ” 7,  and  b^i  = 43.  Then,  ” 33, 
SiQ  =14,  Sg  - 19  and  Sg=19. 


T10 


Figure  2.2.  A Spanning  Tree  for  Example  2 . 1 
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In  the  algorithm,  the  leaving  arc  is  obtained  in  the  following 
manner : 

Let  Q be  the  set  of  nodes  jeVp  with  Sj>0.  A subtree  Tp,  peQ, 
will  be  selected  for  pivot  in  such  a way  that  no  other  node  jeVp  on  the 
path  from  node  p to  node  0 will  have  Sj>0.  Initially,  an  easy  way  of 

determining  Tp  is  to  check  the  depth  functions  of  the  nodes  in  Q and 

select  Tp  such  that  D(p)  = min{D(j):  jeQ) . Throughout  the  algorithm, 
after  a particular  pivot  k,  Tp  is  determined  while  updating  Tk+1.  Let 
Pp'  be  the  set  of  nodes  on  the  unique  path  from  node  p found  in  pivot 
k,  p , to  node  0 in  T^+^-.  Then  new  p is  chosen  to  be  the  node  aePp»n  Q 
with  the  smallest  depth  function.  If  Pp»n  Q = <f> , then  p is  chosen  to  be 

the  node  oceQ  with  the  smallest  depth  function,  and  p'  is  set  to  p.  Once 

a subtree  for  pivoting  is  determined,  the  leaving  arc  is  then  the 
dropped  arc  ((PR(p),p)  in  determining  the  induced  subtree  Tp. 

2.4,5.  Determination  of  the  Entering  Arc 

From  the  construction  of  the  dual  problem,  we  have 

Wi-Wj  < cjj  (i,j)eE.  (2.9) 

Moreover , 

wi-wj  = Cij  if  (i, j)eET.  (2.10) 

Let  Tp  be  the  subtree  selected  for  pivoting.  Also  let 

A - min  { Ajj  - wj -Wi+Cij  ) = Afl  > 0 (f,l)  e E\ET . (2.11) 

ieVp 

j eV\Vp 

Then  we  enter  arc  (f,l)  into  the  tree  T and  drop  arc  (PR(p),p)  out  of 
the  tree  and  update  the  dual  variables  as  follows  : 
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wi  “ wi  + A if  ieVp  (2.12) 

wi  “ wi  if  i€V\Vp . (2.13) 

Lemma — 2.3:  This  dual  solution  is  a basic  feasible  solution  which 
can  be  characterized  by  its  associated  rooted  tree  T*  obtained  from  T 
by  deleting  arc  (PR(p) ,p)  and  entering  arc  (f,l). 

Proof:  We  note  that  for  all  arcs  (z,y)  with  zeVp,  yeVp  or  zeV\Vp, 
y£V\Vp , the  inequality  does  not  change.  Moreover  for  an  arc  (z,y), 
zeVp>  yeV\Vp , the  inequality  holds  due  to  the  selection  of  A.  Finally 
for  the  arcs  (z,y),  zeV\Vp , yeVp  we  have 

wz  - (Wy+A)  < Cij  (2.14) 

which  is  also  satisfied  since  A > 0 . 

Each  time  such  a tree  change  (a  dual  simplex  pivot)  is  performed 
a new  arc  is  added  to  the  tree  which  will  try  to  funnel  the  excess 
suPply  in  Tp  to  a new  node  l<£vp.  Moreover,  since  the  arc  (PR(p)  ,p)  is 
dropped  out  of  T,  it  will  have  zero  flow  in  the  current  attempt  to  find 
the  primal  feasible  solution. 

The  algorithm  to  be  presented  in  the  next  section  will  perform 
the  dual  simplex  pivots  as  follows  : First  a node  p£Q,  with  the 
property  that  for  all  other  nodes  jeVp  on  the  unique  path  from  node  p 
to  node  0 have  Sj<  0,  will  be  selected.  Then  the  arc  (PR(p),p)  will  be 
dropped  out  of  T and  an  arc  (f,l)^E'p  will  enter  the  tree  T.  Here  (f,l) 
is  an  out-of-tree  arc  which  yields  the  value  A calculated  in  (2.11). 

This  process  will  be  repeated  until  - Sj  < 0 for  all  ieVp.  At  that 
point  Q -<f>,  and  the  current  spanning  tree  possesses  an  optimal  solution 
to  the  primal  problem.  The  validity  of  this  algorithm  and  optimality 
results  are  presented  in  Section  2.4.7. 
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2,4,6  The  Dual  Simplex  Algorithm.  DUAL 

In  this  section,  we  present  the  general  algorithm,  DUAL,  for 
solving  the  (HCP) . 

Step  0 . (Initialization')  Construct  the  shortest  path  tree  of  G rooted  at 
node  0 by  using  costs  as  weights.  Let  w^,  ieV  represent  the 
negatives  of  the  shortest  path  lengths  form  node  0 to  node  i. 
Let  T=(V,E-p)  be  this  tree  with  dual  variables,  w^ , ieV. 

Compute  the  depth  function  D(i),  and  by  using  procedure 
PRUNE 1 , ieVp. 

Let  Q - { i: ieVp, S^>0) . 

Set  p’  < 0. 

Set  k < 0. 

Step  1. WHILE  DO 
BEGIN 
k < k+1 

Step  la.  (Determination  of  the  leaving  variable) 

IF  p ' =0  THEN 

Find  Tp  such  that  D(p)  = min  { D ( j ) : jeQ} 

ELSE 

Find  Tp  such  that  D(p)=  min  (D(j):  jeQ  n Pp') 

where  Pp > is  the  set  of  nodes  on  the  unique  path 
from  node  p found  in  iteration  k-1  to  node  0 in  T^. 

END  IF 

Let  (PR(p),p)  be  the  arc  which  is  dropped  to  obtain  Tp. 


Step  lb.  (Determination  of  the  entering  variable) 


Find  A - min  {Ay  - wj -wj+cy  : ieVp,  jeV\Vp,  (i,j)eE)  = Afl 
Let  wj  <---  wj  + A jeVp 
wj  <---  wj  jev\vp 

END 

Step  lc.  (Updating) 

Update  T,Si,  S^.Q. 

IF  Pp  n Q - 4'  THEN 
Set  p’  - 0 

ELSE 

Set  p - p. 

ENDWHILE 

The  current  tree  possesses  an  optimal  solution  to  (MCP) . 

END 

Example  2.2.  The  minimum  cost  flow  problem  with  three  source 
nodes  and  four  demand  nodes  is  solved  by  the  algorithm  DUAL  in  Figure 
2.3. 

20  20 


Figure  2.3  An  Example 
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Initialization 
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Iteration  1.  S4  = 25,  Q - {4},  p=4.  Arc  (0,4)  leaves,  arc  (4,7)  enters, 
A-l. 
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Iteration  2.  Sy=15,  S5-IO.  Q-  (5,7),  p-5 . Arc  (0,5)  leaves,  arc  (4,3) 
enters,  A=l. 
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Figure  2.3  (continued) 
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Iteration  3.  S3  - 10.  Q=-  (3),  p»3.  Arc  (2,3)  leaves,  arc  (4,1)  enters, 
A— 6 . 


10  20  15 


Iteration  4.  Q — <j> . We  have  the  optimum  solution. 


Figure  2.3.  (continued) 


2.4.7  The  Validity  of  the  Algorithm 


The  validity  of  the  algorithm  is  based  on  the  fact  that 
throughout  the  algorithm,  the  flows  on  the  backward  arcs  of  T remain 
positive . 

For  each  node  ieVg,  let  t'i'  = (vj/.E^/)  be  the  subtree  obtained 
from  T by  dropping  the  arc  (i,PR(i))  and  taking  the  subtree  which 
contains  node  i. 

t f 


Let  S 


i “ 2 , , a j 
jevis 


£ b-j 
t * J 

jeviD 


(2.15) 


where  and  are  the  sets  of  supply  and  demand  nodes  in  , 

respectively.  Note  that  for  a node  i€Vg,  is  the  flow  on  the 

backward  arc  (i,PR(i)). 

ft  ft  ft 

Given  an  subtree  = (V^  ,E^  ),  we  define  a subtree 

Ti  - (Vi  ,Ei  ) of  Ti  by  the  procedure  PRUNE2 . This  procedure  is 
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similar  to  procedure  PRUNE1  in  that  we  are  pruning  subtrees  with 
positive  Sj  values.  However  in  PRUNE2 , we  use  the  subtree  T-^ ' as  input 
In  PRUNE2 , given  the  subtree  , we  start  from  the  outermost  nodes  of 
Ti  and  eliminate  the  subtrees  Tj  with  positive  sj  values  and  update 

t t 

the  Sj  values  in  the  remainder  of  the  subtree.  Once  we  eliminate  all 
subtrees  with  positive  Sj  values,  the  remainder  subtree  is  called  Tj_*  ' 
with  its  s[ ' ' value. 


PROCEDURE  " PRUNE2" 

INPUT  : - (Vi',Ei'),ieVB,  D(k),  keVi"  n VF. 

OUTPUT  : A subtree  T-"-  (v£ " , e£ " ) £ T^'  with  the  property  that 
Sj  < 0 for  each  jeVF  n v|"  and  its  S^"  value. 

BEGIN 

INITIALIZE  : Set  T["  <----  T^' , R - ( j € V^'  D VF  : Sj  > 0 ) . 

Zi  - <t>. 

WHILE  R =f=  4>  DO 

Find  jeR  such  that  D(j)  > D(k) , for  all  k+jeR.  Obtain 

Tj“  (Vj-Ej)  • Set  Ti"  “ (Vi"\Vj.  E-"\Ej  U {(PR(j).j)}). 

Zi  <---  Zi  u {j}. 

TT  1 « 9 9 t tit 

Update  Sj  , jeVi  n VB. 

ENDWHILE 

Compute  si"  - S, , , aj  - S bj  (2.16) 

JSVIS  jsV1D 

where  Vig  and  ViB  are  the  sets  of  supply  and 
demand  nodes  in  Ti  respectively. 


END 
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Here  = {j  : jeVp  n , S j >0 } is  the  set  of  nodes  whose 
subtrees  are  pruned  from  subtree  Tj/  with  the  procedure  PRUNE 2 . 

Example  2,3.  Consider  the  minimum  cost  flow  problem  in  Figure 
2.4.  It  is  given  that  ax  - 10,  a2  - 20,  a3  - 20,  a4  = 10 , and  b5  - 5, 
^6  ” 3 , by  — 1,  bg  - 1,  bg  — 0,  b3Q  - 7,  and  b33  = 43.  Then  S2  - 40, 
and  S2  =21. 


Figure  2.4.  A Spanning  Tree  for  Example  2.3 
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At  this  point,  we  will  give  an  overall  review  of  the 
interpretation  of  the  subtrees  T^,  and  T^' ' with  their 

corresponding  S^,  S^,  values  and  their  relationships.  First 

/ t II 

we  note  that  T^,  T^,  S^,  are  defined  for  all  nodes  ieVp,  and  , 

19  9 9 1 III 

T-£  , S , and  are  defined  for  all  nodes  ieVg. 

is  the  subtree  induced  by  the  nodes  ieVp.  Si  value  is  the 
negative  of  the  flow  on  the  forward  arc  (PR(i),i). 

Q is  the  subtree  induced  by  the  nodes  ieVp  when  all  the 

I I I 

subtrees  Tj  Q with  positive  Sj  values  are  pruned  from  T^.  value 
is  the  negative  of  the  flow  on  the  forward  arc  (PR(i).i)  in  Tj_. 

II  II 

is  the  subtree  induced  by  the  nodes  ieVg.  value  is  the 

flow  on  the  backward  arc  (i,PR(i)). 

III  II 

C Tj_  is  the  subtree  induced  by  the  nodes  ieVg  when  all  the 
subtrees  Tj  Q with  positive  Sj  values  are  pruned  from  Tj  . Si" 
value  is  the  flow  on  the  backward  arc  (i,PR(i))  in  Tj ' ' . 

Thus,  we  can  write 

Si  = Si  - 2 S;  for  ieVp  and  (2.17) 

j€Zi 

Si"  - Si'  - 2 S-  for  ieVg.  (2.18) 

jeZi 

where  the  set  Zi  is  the  set  of  nodes  whose  subtrees  are  pruned  from 
subtree  Ti  or  Ti  using  procedures  PRUNE1  or  PRUNE2. 

The  update  of  the  Si,  ieVp,  and  Si  , ieVg,  values  after  a 
particular  pivot  is  nothing  but  flow  update  as  in  the  primal  simplex 
method  around  the  unique  cycle  formed  by  introducing  the  entering  arc 
into  the  spanning  tree.  That  is,  we  add  the  Sp  value  to  all  the  flows  on 
the  arcs  that  have  the  same  orientation  with  the  cycle,  and  subtract  Sp 
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value  from  the  flows  on  the  arcs  that  have  the  opposite  orientation 
with  the  cycle. 

We  will  use  superscripts  ( . )b  and  (.)a  to  indicate  a subtree  (.) 
or  a value  (.)  before  and  after  a particular  pivot,  respectively. 

Let  a pivot  involving  Tp  bring  an  arc  (f,l)^Ex,  feTp,  l^Tp  into 
the  tree.  Let  (x,p)  be  the  leaving  arc  and  let  & be  the  first  common 
node  (join)  between  Pp  and  Ppb.  Consider  the  unique  cycle  formed  by 
introducing  the  arc  (f.l)  into  Tb  (Figure  2.5).  We  can  divide  the  cycle 
into  three  separate  paths  and  examine  the  updating  of  S^,  ieVpb  and 
i€Vgb  values  for  nodes  in  each  path  separately. 


Figure  2.5.  The  Unique  Cycle  Formed  by  Introducing  Arc  (f,l)  into  T 
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For  nodes  iep£\P0b,  if  node  ieVFb , then  ieVFa.  Similarly,  if  node 
i£VBb , then  ieVBa.  That  is  the  orientation  of  the  arcs  on  this 
path  do  not  change  after  the  pivot.  Thus,  we  update  the  values  as 

- Sia  - - Sib  + Spb  ieVFa,  (2.19) 

si’a  - si'b  * Sb  ieVBa.  (2.20) 

Note  that  the  negative  signs  before  SF  values  in  Equation  (2.19) 
is  from  the  fact  that  SF  value  for  a node  ieVF  is  the  negative  of  the 
flow  on  the  arc  (PR(i),i). 

For  nodes  ieP;Lb\P0b,  if  node  ieVFb,  then  ieVFa.  Similarly,  if 
node  ieVgb,  then  ieVga.  Again,  the  orientation  of  the  arcs  on  this  path 
do  not  change  after  the  pivot.  Thus,  we  update  the  values  as 

- Sia  - -Sib  - Sp  ieVFa,  (2.21) 

Si'a  » S-'b  + Sb  ieVBa.  (2.22) 

Note  that  for  nodes  iePfb\Pb,  if  ieVFb  then  ieVBa.  Similarly,  if 
i£vBb>  then  ieVFa.  That  is  the  orientation  of  the  arcs  on  this  path 
changes  after  the  pivot.  Let  v^=  PR(i)  in  Tb.  Then  we  update  the  values 
as 

svi  “ *(Si'b  - Sb)  VieVFa,  (2.23) 

Sv^a  ” -Sib  + Sp  VjeVBa.  (2.24) 

Also  for  the  node  f,  feVBa,  we  have 

sf<a  “ sp-  (2.25) 

Lemma  2,4.  Let  Tp  be  the  subtree  that  is  selected  to  be  pivoted 
on.  Then  Sp  >0.  Moreover  Sp  > S^,  keVp  n VF,  kfp. 

Proof:  From  Equation  (2.17), 

sp  ” sp  + ^ Sk  - E fSk. 

k€Zp  keZp 


(2.26) 
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Since  peQ,  Sp  > 0,  and  2 Sk  > 0,  it  follows  that  Sn  > 0 

keZp  P 

Moreover,  for  any  k+peVp  n VF,  we  have 

^k  ” sk  + 2 Si. 
ieZk 

If  S^.  > 0,  then  since  Zk  u {k}  c Zp , it  follows  that  Sp  > Sk. 
On  the  other  hand,  if  Sk  < 0,  then  Zk  £ Zp  and  thus  Sp  > Sk.  ■ 


(2.27) 


The  following  lemma  states  that  the  flows  on  all  the  backward 
arcs  remain  positive  throughout  the  algorithm. 

~ ma  2'5'  Through°ut  the  algorithm,  the  following  conditions 
remain  true. 

i)  Si  >0  for  all  nodes  ieVg, 

ii)  Si  >0  for  all  nodes  ieVg. 

Proof:  It  is  obvious  that  the  statement  is  true  for  T°  since 

Vb  ‘ * In  T°'  Assu”e  lt  ls  true  at  the  end  of  Iteration  k.  At  iteration 
(k+1)  let  a pivot  involving  Tp  bring  an  arc  (£,1)*ET,  fsTp,  l#Tp  into 
the  tree . 

Since  the  s[  values  of  only  the  nodes  in  the  unique  cycle  shown 
m Figure  2.5  are  changed  by  the  pivot,  we  will  prove  the  conditions 
remain  true  for  only  those  nodes. 

a)  First  we  will  consider  the  nodes  ieCP^b^b  (Figure  2 6) 
From  the  definition  of  the  procedure  PRUNE2 , it  is  obvious  that 
Si"*  - S-"b 

since  Sp  > 0 and  thus  Tp  is  discarded  from  T-" 


in  the  calculation  of 
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Using  Equations  (2.18),  (2.28)  and  the  induction  hypothesis,  we  have 

si"b  - si"a  - Si'a  - 2 Sja  > 0 . (2.29) 

jeZia 

Thus , 

sra>2  s]a>0,  (2.30) 


t t , 

’i 


jeZi£ 


which  completes  the  proof  for  case  a) . 


Figure  2.6  Nodes  between  x and  d before  and  after  the  Pivot 


b)  Consider  the  nodes  ie(Pfb\Ppb)nVpb . (Figure  2.7) 

Letting  Vj_  = PR(i)  , i€(Pfb\Pb)nVpb  in  Tb,  using  Equation  (2.24),  and 
noting  that  v^eVga,  we  have 

va  - -Sib  + Sp 

From  Lemma  2.4,  Sb  > s£,  keVpbnVb,  thus  s"a  > 0 for  v^Vg3. 


(2.31) 
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Before  the  pivot 


After  the  pivot 


Figure  2.7  Nodes  between  f and  p before  and  after  the  Pivot 


To  show  that  condition  ii)  remains  true  on  this  path  after  the 
pivot  consider  a node  ieVpb  and  a node  jeVBb  such  that  i,jePfb\P^  with 
D( j )<D( i) , and  S^a  > 0.  (Figure  2.7)  Let  Vi-PR(i)  and  v«-PR(j) 

j J 

in  Tb . From  Equation  (2.18),  by  dividing  the  set  Zjb  into  two  disjoint 
subsets,  and  by  the  induction  hypothesis,  we  have 


S-"b  - Sj'b  - 2 


k€Zjb\Zib 


2 Skb  > 0 


keZi 


(2.32) 


and  similarly, 


s"  'a 

av. 


s'  'a 

sv. 


sk£ 


kez$  \zva 

i j 


2 Ska 

j 


keZva 


(2.33) 
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We  want  to  show  that  S-^/’3  > 0. 

i 

Note  that  2 sk'a  _ 2 S^b  (2.34) 

keZ^  \Z^a  keZ^b\z[b 

i j J 

since  these  nodes  are  not  in  the  unique  cycle  and  hence  S^a  = S^b. 

From  Equation  (2.24),  we  also  have 

sv'a  “ -Sib  + Spb  . (2.35) 

Subtracting  Equation  (2.32)  from  Equation  (2.33)  and  using  Equations 
(2.34)  and  (2.35),  we  get 


If  A > 0,  then  ' ' a > 0 
i 

hand,  if  A < 0,  then  in 


+ Spb  - 2 S^a  - Sj'b  + 2 S^b  (2.36) 

keZ^a  kez[b 

by  the  induction  hypothesis.  On  the  other 
order  to  show  that  a >0,  it  must  be  true 


that 


_ i t * « t it  K 

sv.  “ sj  b + A > 0,  (2.37) 

l J 

that  is, 

sj'b  - 2 Skb  - 2 S^b  - Sib  + Spb  - 2 S^.a  - S_j'b  + 2 Sfcb  >0  (2.38) 

keZjb\Z-b  kez-b  ke4a  keZ-b 

Simplifying  terms,  we  have  S^/’3  > 0 if 

i 

Spb  - Sib  - 2 Skb  > 2 Ska.  (2.39) 

keZib\z[b  kez^a 
J j 

On  the  contrary,  assume  that 

Spb  * S^b  - 2 Sjjb  < 2 S^3.  (2.40) 

keZjb\Zib  keZ^.a 
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From  our  selection  of  node  j,  we  know  that  Sva  > 0, 

j 


Equations  (2.17)  and  (2.23)  we  have 

£ S^a  - Sa  - - S-'b  + S b. 

ke<a  j 

j 


thus  from 


(2.41) 


Hence 


sj'b  - sib  - s skb  * °-  (2.42) 

keZjb\z[b 

Since  S^b  < £ S^h,  Equation  (2.42)  is  a contradiction  to  Equation 

keZjh 

(2.32).  Moreover,  using  Equations  (2.32)  and  (2.41)  we  have 


Sb  - £ Skb  - £ Skb  > £ S^a . (2.43) 

kez[b  keZib\Zib  keZ^a 

j 

Thus,  we  have  shown  that  Sv  a > 0 if  there  exists  a node  jeVgb  such 

i 

that  i,j,ePfb\Ppb  with  D(j)<D(i),  and  S^a  > 0.  If  S^a  < 0 for  all 

j j 

jGPpb\Ppb  D Vgb,  then  we  can  choose  the  node  j such  that 

D(j)  — min  (D(k)  ; keP^b\Ppb  n Vgb).  Thus  we  have 

£ S^a  - £ S^.a  £ S^b.  (2.44) 

keZ^.a  keZpa  keZp'b\Zjb 

Substituting  Equation  (2.44)  into  (2.43),  we  have 

Sb  - £ S^b  - £ Skb  £ S^b  >0  (2.45) 

kezjh  keZjb\Zib  keZpb\Zjb 


which  agrees  with  Equation  (2.17)  and  the  fact  that  Spb  > 0.  This 
completes  the  proof  for  case  b. 
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c)  In  this  case,  we  will  consider  the  node  feVga.  From  Equation 
(2.25),  and  Lemma  2.4, 

Sf'a  - Sb  > 0.  (2.46) 

To  show  Sf’a  > 0,  we  will  consider  two  subcases  : 

First,  we  will  consider  the  subcase  when  feVp^ . (Figure  2.8)  In  this 
case,  we  can  choose  i— f as  we  did  in  case  b) . Then  from  Equation  (2.43) 
we  have , 

Sb  - Z S^b  2 S^b  - 2 S^a  > 0.  (2.47) 

kez[b  kGZ^b\Zpb  keZ^a 

J j 

Since  Zfa  - Z^a  U Zpb  U Z-jb\Zpb,  and  for  those  nodes  that  are 

j 


Before  the  pivot 


After  the  pivot 


Figure  2.8.  The  nodes  between  f and  p before  and  after  the  Pivot 
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The  second  subcase  is  when  feVgb  (Figure  2.9).  We  select  node  i 

such  that  D(i)-  max  (D(i)  : iG(Pfb\Ppb)nVpb } . 

If  for  each  node  ke(Pfb\P^b)nVBb , Sva  < 0,  where  vk  = PR(k)  in 

k 

Tb,then  choosing  any  je(P£b\Ppb)nVj}b  and  using  Equation  (2.45),  we  have 
Sf"a  > 0. 

However,  if  for  a node  ke(Pfb\Pib)nVgb , S^a  > 0,  then  from  Equations 

k 

(2.17)  and  (2.23),  we  have 


Z Sia  - S5  = -Sk  b + Sb 
, k p 

ieZv3 

k 


(2.48) 


t § t v. 

Since  Sk  u > 0,  we  have 


Sk'b 


Z Si0  > 0 


(2.49) 


iGZkb\Zfb 


i€Zfc 


Substituting  Equation  (2.48)  into  (2.49)  we  have 


Sb  - Z 


Sia 


Z Si 


i«<a  ieZkb\Zfb 
k 


Z S^  > 0. 
ieZfb 


(2.50) 


Since  Sf  a equals  sS,  and  Zfa  = Zkb  U Zva  , the  proof  follows. 
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©... 


> 


Before  the  Pivot 
Figure  2.9.  Nodes  between  f and  p before  and  after  the  Pivot 


d)  Consider  the  nodes  iePFb\P0bn  Vgb  (Figure  2.10).  The  condition 


i)  remains  true  since  from  Equation  (2.22),  we  have 


Si'a  - S['h  + Sp 


(2.51) 


Since  SF  b > 0,  from  the  assumption  and  Sp  > 0 from  Lemma  2.4,  we  have 
S-'a  > 0. 

To  show  that  condition  ii)  remains  true,  we  let  k be  the  node  such  that 
D(k)=  max  { D(t)  , teP^XP^  n VF,  S^b  > 0,  and  s'ta  > 0).  For  all 
those  i on  the  path  P^\P^ , T^  b and  T^  a are  the  same,  therefore 


Sf  a=Si 


< 0 and 


Now,  we  let  j be  the  node  such  that  jePfa\Pkan  Vp,  Sj 


Sja  > 0.  For  all  nodes  w on  Pfa\Pja,  Zwa  - ZWD  U Zfa.  Since 


Ij'  a - Sw  a equals  Sp*5,  it  follows  that 


III  Q t I I U 

c a.  c d 


- s ska  - sf 


> 0. 


35 

( 


(2.52) 


keZfa 

For  nodes  i on  Pja\Pka,  a*S^  ^=*=-Sj':>>0. 


Before  the  pivot  After  the  pivot 

Figure  2.10.  Nodes  between  1 and  9 before  and  after  the  pivot 


What  we  have  shown  with  Lemma  2.5,  is  that  our  selection  rule  of 
the  leaving  arc  maintains  a strongly  feasible  tree,  a feasible  tree 
in  which  every  arc  with  a flow  of  zero  is  a forward  arc,  introduced  by 
Cunningham  [19]  and  Barr,  Glover  and  Klingman  [20].  Moreover,  in  the 
algorithm  we  always  have  positive  flow  on  all  the  backward  arcs. 
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Theorem  2.1:  When  the  algorithm  terminates,  the  final  tree  T 
obtained  by  the  algorithm  possesses  a feasible  (and  optimal)  solution 
to  (MCP) . 

Proof:  Since  for  each  node  ieVg,  we  have  >0  from  lemma  2.5, 
each  backward  arc  in  T*  will  have  a positive  flow.  Each  forward  arc 
that  defines  an  induced  subtree  for  a node  will  have  a non-negative 
flow  from  the  termination  criterion  that  < 0,  ieVp.  Thus,  with  Lemma 
2.1,  the  proof  follows.  ® 

2.4,8  Convergence  and  the  Computational  Bound  of  the  Algorithm 

Definition  2.2:  Given  an  initial  dual  feasible  solution  and  its 
corresponding  tree  T®,  the  initial  total  deficiency  Smax  is  defined  as 
follows : 

^max  “ ^ (2.53) 

ieQ 

During  the  course  of  the  algorithm  Smax  is  defined  in  the  same  manner. 

We  will  also  present  another  way  of  determining  the  Smax  value. 

We  will  use  this  definition  in  the  later  sections.  Noting  that 

Si  + S S-  (2.54) 

jeZi 

if  Si  > 0,  then  the  deficiency  in  a given  subtree  Ti  can  be  expressed 
as  Si . 

Define  a set  K Q Q where  K = { j : Sj  > 0,  and  there  exists  no 
other  node  kePj  with  S^  > 0).  That  is  the  set  K contains  those  nodes 
in  Q which  are  closest  to  the  root  in  every  path.  For  example,  in  the 
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problem  in  Figure  2.11,  if  Q=  { 2, 4, 6, 7, 9),  then  K=  {2,6}.  Then,  we  can 
define  the  Smax  value  as 

Smax  “ 2 Si.  (2.55) 

iGK 


Figure  2.11.  An  Example  for  the  Definition  of  Set  K 

At  iteration  k+1 , let  a pivot  involving  Tp  bring  an  arc  (f,l)^Ei 
into  the  tree.  Consider  the  set  Ni  = {i:  iePifiVp)  . 

Lemma  2,6.  If  all  the  nodes  a e Ni  has  S<^b  < 0,  then  the  Smax 
value  will  reduce  by  at  least  one  after  the  pivot.  On  the  other  hand, 
if  there  exists  a node  « e with  Sab  > 0,  then  Smax  value  will  remain 
the  same  after  the  pivot. 

Proof:  Before  the  pivot,  the  total  deficiency  in  the  subtree  Tp 
is, 

Sp  - 2 Sjh  - 2 ai  - 2 bj 

ieZpb  ievbs  jeVbD 

i 

since  Sp  > 0. 


(2.56) 
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Let  z E such  that  D(z)  - max  {D(°0  } . 

aeNi 

In  proving  Lemma  2.6,  we  will  consider  four  cases. 

Case  a)  Szb  < 0,  but  Sza  > 0. 

The  total  deficiency  in  the  separate  induced  subtrees  Tp  and 
Tz  before  the  pivot  equals 

S Sib  + S Sib  (2.57) 

ieZpb  iezb 

The  total  deficiency  in  the  subtree  Tf  = Tz  U Tp  U { ( f , 1 ) } after 

the  pivot  equals 

ca  ob  i ob 

°z  ~ Dp  T Dz 

Z S[b  + Z S[b  + Szb  (2.58) 

i€Zpb  ieZb 

Since  Szb  < 0,  the  total  deficiency  in  the  combined  subtree  T§ 

decreases  by  an  amount  |Szb|.  However,  since  Szb  < 0,  the  subtree  Tzb 

had  not  been  pruned  from  other  subtrees , T^b  °ceN]_  by  the  procedure 

PRUNE1 . But,  Sza  > 0 indicates  Tza  will  be  pruned  from  other  subtrees 

Taa,  oceN]^.  This  pruning  will  increase  some  S^a  values  by  an  amount  of 

| Szb  | . If  there  exists  a node  with  s4b  > 0,  such  that 

k 

D(ak)  = max  ( D(«)  : oceNi.Sah  > 0},  then  the  contribution  of  node  oc^  to 

the  Smax  value  will  increase  by  the  amount  |Szb|.  Since  T^h  was  pruned 

k 

from  other  induced  subtrees  Ta  before  the  pivot,  and  it  will  be  pruned 
after  the  pivot,  other  Sa  values  will  not  change,  since  Za  = Zb. 


Hence 
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the  total  Smax  value  remains  the  same.  On  the  other  hand,  consider  the 

case  that  for  all  the  nodes  creN^,  Sab  < 0.  Let  al,oc2 “n  the 

subset  of  nodes  in  N]_  with  decreasing  depth  functions  such  that 

Saa  > 0,  i=l n after  the  pruning  procedure. 

i 

We  need  to  show  that  E Saa  < | Sz^  | . 

i 1 


Note  that  Sa 


+ Sr 


n 


n 


n-1 


+ S a 


+ Sg  + 


s'a 


(2.59) 


ieZS 


\z;a 


n 


S<xb  - Sib  + E Sib  + E S[b  (2.60) 

iezb\z|  ieZb 

n 

and 

Sg  - s£  + s£  (from  Equation  (2.21))  (2.61) 

n n ^ 


Since  the  values  for  the  nodes  that  are  not  in  the  unique  cycle  do 
not  change  after  the  pivot,  we  have, 

2 S[b  - E s[a  (2.62) 

ieZb  \Zb  ieZa  \z;a 

n n 

Subtracting  (2.59)  from  (2.60),  and  using  (2.61)  and  (2.62),  we  get 

S S«a  + Szb  - S^b  <0  (2.63) 

1 i n 


indicating  that  the  total  contribution  of  nodes  “l,..-0^  to  the  Smax 
value  is  strictly  less  than  | Sz'->  | . Thus,  we  have  shown  that  for  case  a) 
if  all  the  nodes  cxeN^  have  Sab  < 0,  the  Smax  value  reduces  by  at  least 


one . 
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Case  b)  Szb  > 0 , Sza  > 0. 

Since  Szb  > 0,  zeZzb  and  zeZza.  Thus  from  Equations  (2.57)  and 
(2.58),  we  note  that  the  total  deficiency  in  the  subtree  Tza  remains 
the  same.  Moreover,  Z^^Z^  implies  that  - S^b  for  a=|=z  e N]^.  Thus 
the  Smax  value  remains  the  same. 

Case  c)  Szb  < 0,  Sza  < 0. 

Total  deficiency  in  the  subtrees  Tpb  and  Tzb  equals 

Z S-b  + Sp . 
i6Zzb 

Total  deficiency  in  the  subtree  Tza  equals 
Z s[a  + Z Sia  - Z Sib  + Z s[a 
ieZza\Zfa  ieZfa  ieZzb  ieZfa 

The  equality  in  (2.65)  follows  from  the  fact  that  the  nodes  in 
the  subtree  Tza\Tfa  are  not  included  in  the  unique  cycle. 

From  Equations  (2.25)  and  (2.18),  we  have 

Spb  = Sf'a  - Z S[a  + Sf"a.  (2.66) 

i£Zfa 

From  Lemma  2.5,  we  have  Sf,>a  > 0.  Substituting  (2.66)  into  (2.64)  and 
subtracting  Equation  (2.65)  from  (2.64),  we  see  that  the  deficiency  in 
the  subtrees  involved  in  the  pivot  decreases  by  the  amount  Sf,a. 

I t t « t V 

However  Sf  is  added  to  SaD  value  for  some  oceN]_.  Using  the  same 
argument  as  in  case  a)  if  all  the  nodes  oceNi  have  Sab  < 0,  then  the 


(2.64) 


(2.65) 
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Smax  value  reduces  by  at  least  one.  If  there  exists  a node  oc  with  > 
0 , then  the  Smax  value  remains  the  same . 

Case  d)  Szb  > 0,  Sza  < 0. 

This  case  cannot  occur  because  of  Equations  (2.17),  and  (2.21).  Since 
S2a  = Szb  + Spk,  Szb  value  increases  and  the  Sz  b value  should 
increase  as  well. 

Lemma  2,7.  The  total  number  of  pivots  without  decreasing  the  Smax 
value  is  bounded  by  |V| . 

Proof:  If  the  Smax  value  remains  the  same  at  iteration  k+1,  then 
<xeQk+^  and  D^+ ■*-(«)  < Dk+-*-(f).  Assume  that  <xeQk+l  is  the  node  in 
with  the  lowest  depth  function.  Then,  the  next  pivot  will  be  done  on 
Tak+1  where  Vpk  c Vock+-'-  since  Tak+b  contains  at  least  one  more  node 
than  Tpk.  Thus  the  proof  follows. 

Theorem  2,2.  The  total  number  of  dual  simplex  pivots  of  the 
algorithm  is  bounded  by  M | V | where  M =-  2 a-s  . 

jevs 

Proof:  Since  Smax  < M initially,  and  there  can  be  at  most  |V| 
pivots  before  Smax  is  reduced  by  at  least  one,  the  proof  follows. 

Now,  we  will  state  the  complexity  of  the  algorithm  DUAL  in  terms 
of  arithmetic  operation.  We  note  that  each  pivot  requires  at  most 
0(|V|2)  operations  since  at  most  |V|^  comparisons  are  needed  to  find 
the  value  of  A,  given  by  Equation  (2.11),  at  each  pivot.  However  with 
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the  following  implementation  we  can  reduce  the  complexity  in  terms  of 
operations . 

We  let  kg  be  the  iteration  in  which  Smax  value  has  been 
decreased,  and  kn+^  be  the  next  iteration  in  which  Smax  value  is 
decreased  again.  Thus,  k^,...,kn  are  the  iterations  in  which  Smax 
remains  the  same.  Then  Tpi,...,Tpn  are  the  subtrees  involved  in 
iterations  kp  to  kn.  We  note  that  Tp^  c Tp2  C . . . C Tpn. 

At  iteration  k^ , for  each  node  ieVp^,  we  find 
= min  {Ajj  : jeV\Vp]J  and  store  the  value  of  A^  and  the  node 
jeV\Vpp  which  yielded  the  minimum.  Then  A is  the  minimum  of  the  A^ 
values  for  ieVp^.  Since  we  update  the  dual  variables,  w^  for  iGVp]^,  by 
adding  A,  node  j that  yielded  the  minimum  A^  value  for  each  node  i does 
not  change  for  iteration  k2 , and  A^  values  are  decreased  by  A for 
every  ieVp]^.  Thus,  starting  from  iteration  kp,  if  we  store  the  A^  value 
and  the  name  of  the  node  j which  yields  the  value,  the  only  extra 
comparisons  are  those  to  find  the  A value  which  is  the  minimum  of  the 
Ai  values.  Note  that  we  have  to  replace  a name  and  a value  for  some 
node  i.  However,  each  node  i needs  to  be  "scanned"  only  once  until  the 
Smax  value  reduces  at  iteration  kn+]^.  Thus,  at  most  | V | ^ comparisons 
are  needed  for  iterations  k^  to  kn. 

Theorem  2,3.  The  algorithm  DUAL  requires  at  most  M | V | ^ arithmetic 

operations,  where  M = S a-;  . 

jevs 
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2.5  Scaling  Method  for  Minimum  Cost  Flow  Problems 

In  this  section,  we  scale  the  supplies  and  demands  of  the  minimum 
cost  flow  problem  similar  to  that  of  Edmonds  and  Karp’s  approach  [28], 

We  define  r by  2r‘l  < max  { a^.b^  } < 2r  . (2.67) 

ieVc 

jevD 

In  this  scaling  method,  we  will  solve  r+1  minimum  cost  flow 
problems  each  with  different  updated  supply  and  demand  values.  We  will 
refer  to  each  different  problem  as  a subproblem.  Output  of  one 
subproblem  is  input  to  the  next  subproblem  as  a starting  basic  feasible 
solution. 

In  the  0th  subproblem,  the  supplies  and  demands  are  scaled  as 

aj.0  = fa"i/2^1-  1 and  bj°  «[bj/2*”|  = 1 for  all  ieVs , jeVD  (2.68) 

where  M is  the  smallest  integer  greater  than  or  equal  to  x.  For  the 
zth  subproblem,  the  supplies  and  demands  are  updated  as 

aiZ  = fai/2r '^]  and  bjz  = [bj/2T'z]  for  all  ieVs,  jeVD.  (2.69) 

If  we  use  the  leaving-arc  selection  rule  of  algorithm  DUAL  in 
solving  the  problem  in  a given  subproblem  , it  may  be  possible  to  have 
some  nodes  ieVg  with  < 0 in  the  following  subproblems  because  of 
the  changes  in  a^  and  bj  values.  To  circumvent  this,  the  leaving  arc 
selection  will  be  modified  so  that  Lemma  2.5  holds  true  in  all  the 
subproblems.  That  is,  each  backward  arc  in  the  rooted  tree  T will  have 
a positive  flow  throughout  the  algorithm.  Since  the  algorithm  selects 
only  forward  arcs  as  the  leaving  arc,  with  the  modified  leaving  arc 
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selection,  the  convergence  of  the  algorithm  is  guaranteed.  The  modified 
rule  is  as  follows: 

Assume  that  at  subproblem  z,  Tp  is  the  subtree  selected  to  be 
pivoted  on  by  the  algorithm  DUAL.  Before  determining  the  entering  arc, 
the  scaling  algorithm  first  checks  Sp^,  k > z,  in  the  remaining 
subproblems  with  ap^  and  bj^  as  supplies  and  demands.  If  in  any 
remaining  subproblem  t,  Sp*1  < 0,  then  this  subtree  is  not  selected  for 
pivoting.  (Note  here  that  if  we  pivot  on  Tp  and  the  subtree  remains 
unchanged  until  subproblem  t,  then  we  will  have  a backward  arc  with 
nonpositive  flow  in  subproblem  t.)  Here,  by  remaining  subproblems,  we 
mean  all  future  minimum  cost  flow  problems  with 

ap^  -[ap/2T*kj  and  bj^  -fbj/2T'^]  for  k > z.  (2.70) 

Let  Qz  be  the  set  of  nodes  jeVp  with  S j z , SjZ+l Sjr  > 0. 

Also  let  Kz  C Qz  be  the  set  of  nodes  k such  that  there  exists  no  other 
node  jePk  and  jeQz. 

Lemma  2.7.  Initial  Smax  value  at  the  beginning  of  subproblem  0 is 
less  than  | V | 

Proof:  In  the  subproblem,  we  have  ap1"*  = bj®  = 1 for  ieVg , 

jeVp,  a possibly  unbalanced  minimum  cost  flow  problem,  and  the  initial 
Smax  value  is  bounded  by  |V| . 

Note  here  that  the  unbalanced  problem  does  not  create  any 
difficulties  in  this  algorithm  since  a subtree  Tp  with  Sp  > 0 may  not 
be  pivoted  on  due  to  the  altered  leaving  arc  selection  rule.  Moreover, 
we  may  never  have  a subtree  Tp , with  Sp  > 0,  and  Tp  includes  all  the 
demand  nodes  and  Sp  remains  positive  in  all  the  remaining  subproblems 


45 


since  at  least  in  the  (r)^  subproblem  we  have  a balanced  minimum  cost 
flow  problem. 

Throughout  Sections  2.5  and  2.6,  we  will  use  superscripts  A and  B 
in  the  following  manner.  (.Z)A  indicates  the  subtree  (.)  or  a value  (.) 
at  the  end  of  subproblem  z,  and  (,Z)B  indicates  the  subtree  (.)  or  a 
value  (.)  at  the  beginning  of  subproblem  z. 


For  any  subproblem  z > 1,  we  have 


2aiz"1  - 1 < aiz  < 2aiz"1 

ievs 

(2.71) 

2bjz*1  - 1 < bjz  < 2bjZ'1 

jeVD. 

(2.72) 

Note  that  (S§)B,  the  Sa  value  oceVp  at  the  beginning  of  subproblem  z, 
equals  twice  (S^2’^ )A, the  Sa  value  at  the  end  of  subproblem  (z-1),  if 
all  bjz  and  a^z,  ieTag , jGTaD  assume  their  largest  possible  values  in 
(2.71)  and  (2.72).  In  general,  we  can  write 

(SZ)B  - 2(S4Z-1>)A  + |(ND)Z|  - |(NS)Z|  (2.73) 

where  (ND)g  - { j : jeVaD  and  b|  - 2bjz'1  - 1 ) and  (Ng)z  = { i : i eVaS  and 
ap2  - 2aiz'1  - 1 ) . 

With  the  following  lemma,  we  will  show  that  the  S^ax  value,  the 
total  deficiency  at  the  beginning  of  subproblem  z is  bounded  by  a 
polynomial  of  |V| . However,  before  stating  the  lemma,  we  will  consider 
the  case  when  for  a node  oceVp,  (S«Z‘1)A  > 0.  Since  (SaZ‘-*-)A  > 0,  it 
implies  that  there  exists  a subproblem  t > z-1  where  (ScJct)B  < 0. 
Otherwise,  node  oc  would  be  in  Qz*^  and  subproblem  (z-1)  would  not  be 
terminated.  Therefore,  even  though  (Saz)B  may  be  positive  at  the 
beginning  of  subproblem  z,  we  will  not  pivot  on  that  subtree,  since 
Qz,  thus  (Saz)B  will  not  be  added  to  the  S^ax  value  for  that 
subproblem.  In  such  a case,  in  order  to  eliminate  positive  (S^Z)B  in 
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calculating  the  S^ax  value , we  can  temporarily  add  a dummy  demand  node , 

d,  to  T«z  with  b^=  (S^2)®,  thus  making  T„z  not  eligible  for  pivot  at 

the  beginning  of  subproblem  z . After  adding  the  dummy  demand  node , we 

can  recalculate  the  (S^2)®  values  for  those  nodes  keVp  on  the  path  P,*. 

At  the  end  of  subproblem  z,  we  can  remove  the  dummy  node  d.  Keeping 

this  discussion  in  mind,  we  can  define  S^ax  value,  the  total  deficiency 

at  the  beginning  of  subproblem  z as  follows: 

Smax  “ 2 (s[z)B  (2.74) 

ieQz 

or 

Smax  - 2 [ (S±Z)B  - 2 [(sjz)  : jeTi(  S^z  > 0 but  j<j=  Qz]]  (2.75) 

ieKz 

Lemma  2.8.  Sj^ax,  total  deficiency  at  the  beginning  of  subproblem 
z,  z > 1 , is  bounded  by  |Vj)j  . 

Proof : In  this  proof,  we  will  consider  three  cases. 

If  aeKz,  and  (S§*1)A  < 0,  then  from  Equation  (2.73),  (SZ)B 
increases  the  most  if  (Nj))z  = VaD,  and  (Ng)z  = <f> . Thus  (S§)B  is  bounded 
by  |VaD| . 

If  Kz,  and  (S«*^)^  > 0,  but  (Saz‘l)A  < 0,  from  Equation  (2.17) 
note  that 

(Sa‘1)A  “ (S«Z*1)A  + 2 (SiZ‘1)A.  (2.76) 

ie(Zz‘1)A 

Substituting  Equation  (2.76)  into  (2.73),  we  have 

(SZ)B  =2(S;z'1)A  + 2 2 (S-Z'1)A  + |(ND)Z|  -|(NS)g| 


(2.77) 
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Even  though  the  value  2 E (S^Z‘^)A  may  be  positive,  it  will  not  be 

ie(z2-1)A 

added  to  the  S^ax  value.  Thus  (S§)B  value  that  will  contribute  to  the 
smax  value  increases  the  most  if  (ND)g  - Vo^,  and  (Ng)g  - <f> . Thus  (Sg)B 
is  bounded  by  |Vaj)|  . 

If  for  a node  cc,  (S„"1)A  > 0,  and  (Sa2'^)A  > 0,  it  implies  that 
there  exists  a subproblem  t such  that  < 0,  thus  a ^ Qz  or  Kz. 

Since  for  i,jeKz,  V^nVj  - and  S^ax  is  the  sum  of  (SZ)B  values 
for  ocgK2  it  follows  that  Sjjjax  < |Vq|  . • 

Corollary  2.1.  For  ieQz , (S^Z)B  is  bounded  by  |V^g|. 

Proof : From  Lemma  2.8,  we  know  that  for  a node  <xeKz , (SZ)B  is 
bounded  by  |Va£)|.  The  subtree  (T^)B  is  the  union  of  the  subtrees 
(T^z)B,  i€QznTz.  Since  ieQz , we  know  that  (s[z"-'-)A  < 0.  Thus,  those 
demand  nodes  in  (ND)gn(T[z)B  that  increase  (Sg)B  value  will  also 
increase  the  (S^Z)B  value  as  well.  ® 

Theorem  2.4.  Total  number  of  pivots  in  the  scaling  algorithm  is 
bounded  by  (r+l)|V|^.  Moreover,  the  algorithm  requires  0[(r+l)|V|B] 
arithmetic  operations. 

Proof:  Since  we  solve  (r+1)  subproblems  each  with  initial 
deficiency  of  at  most  |V|  by  using  Theorems  2.2  and  2.3,  the  proof 
follows.  * 
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2.6  A Different  Scaling  Algorithm  for  Minimum  Cost  Flow  Problems 

In  the  scaling  algorithm  presented  in  Section  2.5,  for  each  pivot 
at  subproblem  z we  have  to  calculate  the  S^,  ieVp  values  for  all  the 

remaining  subproblems  z+1 r , and  this  increases  the  computational 

time  for  each  pivot. 

In  this  section,  we  present  another  scaling  algorithm  in  which 
the  leaving- arc  selection  rule  is  changed  so  that  we  would  require  less 
computational  effort  at  each  pivot.  In  the  new  arc  selection  rule  only 
the  S-£  values,  ieVp,  for  the  current  subproblem  z and  the  Tth 
subproblem  are  considered.  Since  we  do  not  consider  the  intermediate 

subproblems  (z+1 r-1) , we  may  have  some  backward  arcs  with  a 

nonpositive  flow  in  an  intermediate  subproblem.  That  is,  if  we  pivot  on 
a subtree  Tp  in  subproblem  z while  Sp,r  > 0,  but  S^  < 0 for  subproblem 
t,  z < t < r,  and  if  Tp  remains  unchanged  until  subproblem  t,  we  may 

t I I *- 

have  Sp  < 0 in  subproblem  t,  which  is  a contradiction  to  Lemma  2.5 
which  states  that  all  S^  >0  throughout  the  algorithm.  Thus,  for  any 
subproblem  z=fr , we  can  have  subtrees  T^ , ieVp,  that  are  not  strongly 
feasible.  The  only  subproblem  in  which  the  subtrees  remain  strongly 
feasible  is  the  rth  subproblem. 

In  a strongly  feasible  tree,  the  deficiency  caused  by  a node  ieVp 
was  simply  its  Sj_  value  when  S^>0.  However,  when  there  exists  some 
backward  arcs  with  nonpositive  flows  in  the  tree,  then  the  deficiency 
caused  by  a node  i has  to  be  redefined.  In  order  to  define  this 
deficiency  at  subproblem  z,  we  will  define  a subtree  T^2  and  its 
associated  S^z  value  for  ieVp. 
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Let  Qpz  be  the  set  of  nodes  with  S^z  > 0,  ieVp. 

Given  a subtree  T^z  - (V£Z,E^Z)  for  ieQpz , we  define  a 
subtree  Tp2-  (V£Z,E^Z)  of  T^z  by  the  following  procedure,  PRUNE3 . In 
PRUNE3,  starting  from  the  tip  nodes  of  subtree  T^z,  we  eliminate  the 
subtrees  Tj  z,  jeV-^2  n Vg  with  Sj  2 < 0. 

PROCEDURE  " PRUNE3" 

INPUT  : Tiz  - (V£Z , E^z) , for  ieQp2 , and  D(k) , kev[z  n VB. 

OUTPUT  : A subtree  T^z-  (Vj_z,Ej_z)  Q T^z  with  the  property  that  for 
each  jeViB,  Sj  z > 0,  and  its  associated  S^z  value. 

BEGIN 

INITIALIZE  : Set  Tiz  < T^z , R = { j e v[z  n VB  : s]  "z  < 0 ] . 

Z^  = 0 . 

WHILE  R + 0 DO 

Find  jeR  such  that  D(j)  > D(k)  , for  all  k=j=jeR. 

Obtain  Tj  "z=  (Vj  ' 'Z,E- ' ' z)  . 

Set  Tj2  - (ViZ\Vj"z,  EiZ\Ej"z  U { (j  , PR(J  ) ) } ). 

Let  M C Z^z  be  the  set  of  nodes  m in  Zp2  whose  unique 
path  to  node  0 passes  through  node  j . 

Set  Ziz  <---  Ziz  U {j } \ M 
Update  Sj’’2,  jeV^n  Vg  and  R. 

ENDWHILE 

Compute  Siz  = S aj  E bj  (2.78) 

j^iS  Je^iD 

where  V^g  and  V^g  are  the  sets  of  supply  and  demand 
nodes  in  Tj_ , respectively. 


END 
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Note  that  the  subtree  T^z  C T^z  is  obtained  when  all  the  backward 
arcs  with  nonpositive  flows  are  removed  from  Tpz . Z^z  is  the  set  of 
nodes  j , jeVg  whose  subtrees  are  pruned  from  subtree  T^z  in  subproblem 
z . Thus , 

Siz  = S[z  + E |Sj"z|  for  ieQlZ-  (2.79) 

jeZiZ 

If  for  a node  ieQpz , there  exists  some  nodes  jeVpn  Vg  such  that 

I t t ~ 

Sj  < 0,  then  we  define  the  contribution  of  node  ieQp2  to  the  total 
deficiency  to  be 

^max.i  = siZ  +\^iZ\-  (2.80) 

On  the  other  hand,  if  for  a node  ieQ^2  there  exists  no  node  jeV^n  Vg 
such  that  Sj  2 < 0,  then  the  contribution  of  node  i to  the  total 
deficiency  is  S2axi  - Sp2  - s[z. 

One  may  interpret  the  S^ax^  value  for  ieQp2  as  the  deficiency  of 
an  analogous  strongly  feasible  subtree  obtained  from  subtree  Tp2  by 
replacing  each  subtree  Tj  z,  jeZ-^,  with  a supply  node  that  has  a 
weight  of  1. 

When  we  remove  all  subtrees  T^z  for  ieQp2  from  the  tree  Tz , we 
have  a remaining  subtree  that  has  no  negative  flows  on  forward  arcs , 
but  may  have  some  nonpositive  flows  on  backward  arcs.  We  will  call  this 
remaining  subtree  Tp.  There  may  be  some  nodes  ieTpnVp  such  that 
^max.i  ” Siz  + E | S j z|  + |Z£Z| 

jeZiz 

can  be  positive.  Here  Z^2  is  the  set  of  nodes  in  T^z  with  S j ' ' z < 0. 

Let  Q2Z  be  the  set  of  nodes  ieTp  n Vp  satisfying  this  condition. 
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We  can  find  the  nodes  in  Q2Z  in  the  following  manner  : 

Set  Q2Z  — 0. 

Let  R = { j e TR  n Vg  : Sj ’ ' z < 0 } . 

WHILE  R + 0 DO 

Find  jeR  such  that  D(j)  > D(k)  for  all  k=fjeR. 

Find  the  node  iePjnVp  such  that  |S^Z|=  min{|SR'z|,  iePjfWp)  . 
Use  Procedure  PRUNE3  as  T^Z  as  the  input  to  find  T^2 . 

If  Siz  - S Sj  2 + |ZiZ  | > 0 , remove  Tj  from  TR  and 
set  Q2Z  - Q2Z  U { i } . 

Update  TR,  R,  S^2  for  ieTRnVp,  and  Sj’’2  for  jeTRnVB. 
ENDWHILE. 


Again  we  can  interpret  the  S^axi  value  for  ieQ2Z  as  the  deficiency  of 
an  analogous  strongly  feasible  subtree  obtained  from  subtree  T^Z  by 
replacing  each  subtree  T j ' ' 2 , jeZj^  with  a supply  node  that  has  a 
weight  of  1. 

If  we  let  QZ  = Qp2  u Q2Z , then  we  can  define  the  total  deficiency 
in  subproblem  z , S^axz  as 

^max  Z “ 2 S^ax  ^ (2.81) 

i€Q2 

Example  2,4.  Consider  the  problem  in  Figure  2.12  with  5 source 
nodes  and  8 demand  nodes.  It  is  given  that  ap  = 4,  a2  = 5 , 33  = 4, 

34  = 11,  and  35  - 10,  and  bg  = 5,  by  = 4,  bg  = 2 , bg  = 3 , bpo  = 5, 
bn  = 3,  bp2  - 1,  and  bp3  = 11.  Here  Sp2  - S]^'  - 1,  Sp2  “ 4. 

Z12  = (2,3)  and  SraaX)p2  — 6. 
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max , i 


Figure  2.12.  Example  for  the  Definition  of  S 
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2,6.1.  The  Scaled  Dual  Simplex  Algorithm.  DUALSC 

(Initialization)  Calculate  r by  2r‘^  < max  { a^,bi  ) < 2r  . 

ieVs 

jeVD 

Construct  the  shortest  path  tree  of  G rooted  at  node  0 by  using 
costs  as  weights.  Let  w^,  iGV  represent  the  negatives  of  the 
shortest  path  lengths  form  node  0 to  node  i.  Let  T=(V,E'p)  be 
this  tree  with  dual  variables,  w^,  ieV. 

Set  z < 0. 

Set  k < 0. 

WHILE  z < r DO 
Calculate 

a^2  =[a^/2r_2l  and  bj2  = fbj/2T_z]  for  all  ieVg , jeVj). 

For  each  node  ieVp,  compute  D(i),  S^7"  values  and  S^ax  ^ 
values 

Let  Qz  = { i : ieVp, S^aXi ^ > 0 and  S^r  >0). 

Set  p < — 0. 

WHILE  Qz+  o DO 
BEGIN 

k < ---  k+1 

(Determination  of  the  leaving  variable) 

IF  p'=0  THEN 

Find  Tp  such  that  D(p)  = min  { D ( j ) } 

jeQz 


ELSE 
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Find  Tp  such  that  D(p)=  min  { D ( j ) } 

jeQznPp' 

where  Pp  is  the  set  of  nodes  on  the  unique  path 
from  node  p found  in  iteration  k-1  to  node  0 
in  T^. 

END  IF 


Let  (PR(p),p)  be  the  arc  which  is  dropped  to  obtain 


(Determination  of  the  entering  variable) 


Find  A — min  {Ajj  — wj -w^+c-jj  ) = Af]_ 
ieV„ 


(i.j)el 

) 

Let  wj 

< — 

Wj  + A 

jeVp 

Wj 

< — 

Wj 

jeV\Vp 

END 


(Updating) 


Update  T, 

^max , 

• z s 
1 > b 

IF  Pp  n Qz 

= 0 

> 

THEN 

Set  p 

= 0 

ELSE 

Set  p 

" p. 

END  IF 


ENDWHILE 

z < z+1 

ENDWHILE 

The  current  tree  possesses  an  optimal  solution  to  (MCP) . 


END 
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2.6,2  Computational  Bound  of  the  Algorithm 

In  the  following  lemma,  we  assume  that  we  are  at  subproblem  z, 
and  we  have  done  k pivots  at  subproblem  z.  To  eliminate  confusion 
between  pivot  indices  and  subproblem  indices,  we  omit  the  subproblem 
indices.  Superscripts  a and  b are  used  to  indicate  a value  after  and 
before  the  pivot,  respectively. 

At  iteration  k+1,  let  a pivot  involving  Tp  bring  an  arc 
(f,l)  ^ Ep  into  the  tree.  Consider  the  set  Np  = { i:  iePpnVp)  in  . 

Lemma  2.11.  If  all  the  nodes  in  oceNp  have  Smax  oc^  < 0,  then  the 

I 

smax  value  decreases  after  the  pivot.  On  the  other  hand,  if  there 
exists  a node  « £ Nj  with  Smax<xk  > 0,  then  the  Smax  value  may  remain 
the  same  after  the  pivot. 

Proof:  One  way  of  proving  the  lemma  is  to  use  the  analogous 
strongly  feasible  tree  obtained  by  replacing  the  subtrees  with  negative 

til 

Sj  values  by  supply  nodes  of  weights  one.  Since  we  can  interpret  the 
Smax  value  of  the  original  tree  as  the  Smax  value  of  this  analogous 

strongly  feasible  tree,  we  can  use  the  same  proof  as  in  Lemma  2.6. 

However,  we  include  the  proof  of  the  lemma,  using  the  original  tree. 

In  this  proof,  we  will  first  consider  the  change  in  the  Smax,i 
values  of  nodes  i in  the  subtree  Tp^.  Afterwards,  we  will  consider  the 
change  in  the  Smaxp  values  of  the  nodes  i in  N]_.  For  all  those  nodes 
ieVp  that  are  not  on  the  unique  cycle,  the  Smaxp  values  will  not 
change . 

In  the  subtree  Tp^,  only  the  Smaxp  values  of  the  nodes  i that 
are  on  the  path  Pf^\Pp^  will  change  by  the  pivot  since  these  are  the 
only  nodes  that  lie  on  the  unique  cycle. 
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Let  v1(...,vK  be  the  nodes  in  Pfb\PpbnQb  such  that  there  exists 
a node  j € n (Pfb\Ppb) . (Figure  2.14) 

These  are  the  nodes  in  Qz  whose  induced  subtrees  contain  a node  jeVg  on 
Pfk\Ppk  with  Sj  b < 0.  We  will  call  this  node  j (vj) . The  case  when 
there  does  not  exist  such  nodes  will  be  handled  subsequently. 

Let  ti  be  the  last  node  of  Tj (v£)  " 'b  on  the  path  Pfb\Ppb.  That  is, 
D(ti)=max{D(t)  : tePfb\Ppb  n VjCv^h. 

We  will  define  a subtree  T^b  as 

rqb  - U Tib. 

ie(Pvb\Ppb)oQb 

In  general  rkb  is  defined  as 


rkb  = U 

ie(P  b\Ptb  )r*)b 

k k-1 


Tib  for  k=2 K. 


Note  that  the  subtrees  rkb  and  rk+ib  are  separated  by  the  subtree 

Tj  (vk)b  • Tbe  contribution  of  the  nodes  in  Tkb  to  the  value  can  be 

expressed  as 

2 Sib  + z |Zib| 

ierkbnQb  ierkbnQb 


— 2 

ierkb  nvs 


2 bj  + 2 

jerkbnvD  ierkbnQb 


I- 


We  will  let  qk 


PR(j(vk))  for  k-1 K in  Tb . 
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Before  the  Pivot 


After  the  Pivot 


Figure  2.13.  Nodes  in  Pf^\Pp^  before  and  after  the  Pivot 
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After  the  pivot,  we  define  T^a  as  the  subtree  rooted  at  node 
such  that 


= U T^3 

ie(Ppa\Pqia)nQa 

The  contribution  of  the  nodes  in  T^a  to  the  S^,axa  value  can  be 
expressed  as 

E Sia  + S |Zial 

i€rianQa  ierianQa 


‘ E a^  * S bj 

ieri^Vs  jer^nvo 


E | Zia 
ieriaoQa 


Since  the  node  sets  of  Tpa  and  are  the  same,  we  have 


2 - 2 bj 

ier^nVg  jer!anVD 


ier;Lbnvs 


E bj 

jer^nvo 


Moreover , 

2 I Zfa  I - E | Zib  | - 1 

ier1anQa  ierj^hnQh 

since  the  node  j (v]_)  is  no  longer  in  Z^a.  Thus  the  contribution 
of  the  nodes  in  T^a  to  Smax  value  decreases  by  one  after  the  pivot. 
We  will  now  look  at  the  nodes  in  Tj (v  )b  after  the  pivot. 

Since  the  node  sets  of  Tj  (v^)b  anc*  Tt  a are  the  same,  we  have 

Sj  (▼,)' 


< 0. 
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Thus  subtree  Tt^  3 will  be  pruned  by  procedure  PRUNE3  after  the  pivot. 


In  general,  we  can  define  the  subtree  rk3  k=2,...,K  as  the 
subtree  rooted  at  node  qk  as 


rk3  - u T^3. 

iePt  a\PaanQa 

k-1  qk 

The  contribution  of  the  nodes  in  I\a  k=fl  to  the  S^ax3  value  is 

expressed  as 

E Si3  + S | Z±a | 

ierkanQa  ierkanQ3 

= s aj_  - 2 bj  + 2 I Z±a  I . 

ierk3nvs  jerkanvD  ierkanQ3 


Since  the  node  sets  of  Tk3  and  Tkb  are  the  same,  and 


2 | Zi3  | 

ierkanQ3 


ierkbnQb 


I Zlb  i . 


the  contribution  of  nodes  in  rk3  to  S^ax3  value  equals  to  the 

contribution  of  nodes  in  Tkb  to  S^axb  value. 

The  value  Sf  3 can  take  two  possible  values  depending  on  whether 

node  f was  in  the  subtree  Tkb  or  in  the  subtree  Tj  (y  )b  in  Tb . The 

K 

change  in  the  contribution  of  nodes  in  Tpb  to  the  S^ax  value  is 
reflected  in  the  Sf  3 value.  We  will  consider  each  case  separately. 
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Case  1.  If  f was  in  the  subtree  Tj )b, 

1C 


sf 


Sj(vK)b  < 0. 


then 


In  this  case  the  total  contribution  of  the  nodes  in  Tpb  to  the  S^ax 
value  decreases  by  one  since  for  each  pair  of  subtrees  T^3  and  , the 
contribution  of  nodes  in  the  subtrees  remains  the  same  for  k=2,...,K, 
and  reduces  by  one  for  k=l . 

Case  2.  If  f was  in  the  subtree  then  we  have 

S V + Sf"3  - 2 aL  - S bj  - 2 Sib. 

ierKanQa  ierKanvs  jerKanvD  i<=rKbnQb 


* t t « 

Thus  Sf  is  the  difference  between  the  contribution  of  the 
nodes  in  T^b  and  T^3.  Here  we  want  to  note  that  if  we  cannot  find  any 
nodes  Vf,  then  we  will  have  one  subtree  Tfb  and  Tfa  with  K=1 . 

Now  we  will  look  at  the  change  in  the  contribution  to  the  S^x 
value  of  the  nodes  on  Nf.  According  to  the  sign  of  Sf’3,  we  will  again 
consider  two  cases. 

Case  1 . Sf ' ' 3 < 0 . 

Note  that  Sf  3 < 0 implies  that  the  contribution  of  nodes  in  Tpb 


t0  ^max  value  has  already  been  reduced  by  one . 

Case  la)  If  there  exists  a node  qeNf  such  that  SmaXqb  > 0,  and 
D(q)“  max{D(j ) : Smaxjb  0},  then  Sf"3  value  will  be  added  to  Sq'b 
value  and  Sj  b values  for  j6Tq  b n Pf,  decreasing  both  values  by 
the  same  amount.  Using  Equation  (2.79)  we  can  see  that  Sq  value 
remains  the  same.  However, 

lZq  I lzqbl  + since  Zq3  - Zqb  U (f).  Thus,  Smaxqa  = Smaxqb  + 1. 
Since  the  contribution  of  nodes  in  Tpb  to  Smax  value  is  reduced  by  one 
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and  the  contribution  of  node  q to  Smax’  value  is  increased  by  one, 
the  Smax  value  remains  the  same. 


Case  lb)  On  the  other  hand,  if  all  the  nodes  qeNf  have 
^max , q ^ then  Sjjjax^q  “ ^max,qb  1 — 0*  Thus,  node  q has  still  no 
contribution  to  Smax  value.  Since  the  contribution  of  nodes  in  Tpb  to 

I t 

Smax  value  is  reduced  by  one,  the  Smax  value  is  reduced  by  one. 


Case  2.  Sf  a > 0. 

Note  that  Sf  a > 0 implies  that  the  contribution  of  nodes  in  Tpb 

t III 

to  Smax  value  has  already  been  reduced  by  Sf  a > 0. 

Case  2a)  In  this  case,  again  if  there  exists  a node  qeNf  such 

that  smax,qb  ^ °.  and  D(q)-  max{D(j)  : Smaxjb  > 0),  then  the  Sf"a 

value  will  added  to  the  Sq'b  value  and  Sj ' 'b  values  for  jeTq'bnPfnVB, 
increasing  both  values  by  the  same  amount,  jf  there  exists  a node 
jeTq  bnPfnVg  with  Sj  b < 0,  then  we  can  have  two  subcases: 

i)  If  Sf  a > | S j b| , then  the  Sj  a value  will  be  positive,  and 
the  Sq  value  will  increase  by  the  amount  Sf’’3  - | S j * ' b | . This  follows 
from  the  fact  that  in  Equation  (2.79),  Sq’  increases  by  Sf’,a, 
and  j<£Zqa.  Since  the  contribution  of  nodes  in  Tp  to  the  S^ax  value 

has  reduced  by  the  Sf  a value,  the  combined  contribution  of  the  nodes 

in  Tp  and  node  q to  the  Smax  value  will  decrease  by  the  amount  of 

isrbi- 


ii)  If  Sf  a < | S j b|,  then  the 
in  Equation  (2.79),  the  Sq  and  Sj ’ ’ 
while  j remains  to  be  in  Zq.  Thus, 


Sq  value  remains  the  same,  since 
values  increase  by  the  same  amount 
the  Smaxq  value  remains  the  same. 
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Since  the  contribution  of  nodes  in  Tp  to  the  Smax  value  has  been 
reduced  by  the  Sf  a value,  the  combined  contribution  of  the  nodes  in 
Tp  and  node  q to  the  Smax  value  will  decrease  by  the  amount  of  Sf  a. 

If  there  exists  no  node  jeTq  bnP]_nVg  with  Sj  b < 0,  then  the 
Smax.q  value  increases  by  the  amount  of  the  Sf  a value,  thus  the 
Smax  value  remains  the  same . 


Case  2b)  If  there  does  not  exist  a node  qeNf  such  that 
£>max,qk  > 0,  however  there  exists  some  nodes  qeNf  such  that 
Smax,qa  - 0 , we  will  consider  the  following: 

Smax,qk  < 0 implies  that,  Sq  b < 0,  and  Z qb  = 0 for  all  qeNf. 

This  means  that  there  are  no  backward  arcs  with  nonpositive  flows  on 
the  path  Pf . 

, Let  t eNf  such  that  D(t)  = max  (D(q) : Smax  qa  > 0). 

t t I « — 

Since  Sf  d is  added  to  all  the  flows  on  backward  arcs,  Zta  = 0.  Thus, 
Smax,ta  - Sta  = Stb  + Sf  a < Sf  a. 

Since  the  contribution  of  node  t to  the  Smax  value  increases  by  the 
value  Stb  + Sf  a , where  Stb  < 0,  the  combined  contribution  of  the 
nodes  in  subtrees  Tp  and  Tt  to  the  Sraax  value  decreases  by  the  amount 
of  |Stb| . Using  the  same  argument  as  in  Lemma  2.6  case  a),  we  can  show 
that  the  increase  in  Smax  oc  values  for  <xePt  is  strictly  less  than 
|Stb| . Thus,  the  Smax  value  decreases  by  at  least  one.  This  completes 


the  proof  of  the  lemma. 
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Lemma  2.12.  The  total  number  of  pivots  without  decreasing  Smax 
value  is  bounded  by  |V| . 

Proof : Same  as  Lemma  2.7. 

Note  that  in  this  new  scaling  algorithm,  DUALSC,  we  do  not  pivot 
on  the  backward  arcs  with  negative  flows  since  in  the  original  problem 
we  know  that  there  is  a positive  flow  on  that  arc. 

To  find  the  computational  bound  of  this  new  scaling  algorithm  we 
consider  the  following  lemmas.  In  Lemma  2.13,  we  will  give  a lower 

t t I — 

bound  on  the  Sj  value  for  a node  jeVg. 

Lemma  2.13:  At  the  beginning  of  subproblem  z,  the  Sj  2 values 
for  all  jeVg  is  bounded  from  below  by  -|(Vjg  )z|  . 

Proof:  If  for  a node  jeVg,  Sj  z < 0,  then  we  will  show  that 


> ' I V j D ' 2 I • 

If  Tj  z is  the  same  subtree  as  Tj  z+ ^ Tj  r,  that  is 


Z j z=Z j z+l= . . . =Z j r even  though  we  change  the  supplies  and  demands  in 


each  subproblem  (Figure  2.14),  then  we  can  write 


(Sj"2+1)  - 2 (Sj  " z)  + |(NjD)z+1|  - |(NjS)z+1| 


(2.82) 


where  (Njg)z+-'-  = { k:keVjg  z+ and  b^2*^  = 2b]^z  - 1 } and 
(NjS)z+1  - { i:ieVjs'z+1  and  aiz+1  = 2aiz  - 1 ). 

In  such  a case,  the  minimum  value  (Sj  z)  can  assume  is  bounded  by 
- |VjD  |.  To  see  this,  assume  on  the  contrary  that 
(Sj  z)  < -|Vjg  z|.  Then  from  Equation  (2.82),  the  most  increase  in 

t I t T.  lit _ 

Sj  , k=z r,  one  can  achieve  is  when  (Njj))z  = Vjg  2 and 


(Njs)Z  = 0-  Thus 


(Sj"2+1)  - 2(Sj  * ' z)  + | Vjo' z | < -|vjo'Zl, 


(2.83) 
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and 


<srZ+2>  - 2<si"2+1)  + |vjD'Zl  ^ • 

Continuing  in  this  manner,  we  get 

(Sj  T)  < -|VjQ  z|,  contradicting  to  the  fact  that 

t t t _ 

Sj  >0,  for  all  jeVg,  that  is  in  the  original  problem,  we  have 
negative  flows  on  the  backward  arcs. 


(2.84) 


no 


However,  when  zj  z+Zj  z+1=f  ■ ■ -+Zj  T , then  the  subtree  Tj  ' ' changes 
during  the  remaining  subproblems.  Keeping  this  in  mind,  we  will 
consider  four  cases  shown  in  Figure  2.15. 

Case  1 : Consider  a node  ieTj ’ such  that  ieZjz  and  ieZj T . (For  example 
node  i in  Figure  2.15a). 

Case  2 : Consider  a node  ieTj ’ such  that  ieZjz  but  i^Z j r . 
node  i in  Figure  2.15b). 


(For  example 
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We  will  use  the  same  proof  for  cases  1 and  2. 

On  the  contrary,  we  will  assume  that,  Sj  z < - | Vj j)  z|.  We  will 

consider  the  following  subproblems,  t]_ , , t2  , q2  > • ■ • tK> 

where  z < t^  < q^  C.  . . < qj<;  < r . In  subproblems  z to  t^-1,  we  have 

ieZj . In  subproblems  t^  to  q^-1,  we  have  i^Zj . In  general,  in 

subproblems  t^  to  q^-l,  we  have  i^Zj , and  in  subproblems  q^  to 
(tk+i  -1),  we  have  ieZj . 

Since  the  subtree  does  not  change  until  subproblem  t]_,  as  in  the 
previous  case  we  can  show  that  at  subproblem  tj_-l, 


(S; 


V1) 


< - V 


jD 


(2.85) 


At  subproblem  tj_,  subtree  Tj.^]_  is  included  in  Tj  t]_ . Thus  we  know  that 
< 0,  and  the  most  increase  in  (Sj  *-]_)  is  achieved  by 


ti  . ,c:"Vi  + | v j d ' z I + siV 


(2.86) 


In  the  remaining  subproblems  until  subproblem  q]_,  where  subtree  T^ 
remains  to  be  a part  of  subtree  Tj  , nodes  in  Iv^q  | help  in 
increasing  Sj  value  and  value  if  they  assume  their  lower  bounds. 
Thus,  the  most  increase  in  Sj  ti+'*'  is  achieved  when 


V1 


2S- 


#t. 


|VjD 


tl+li 


(2.87) 


Realizing  that  Vj  q' = Vjj)’z  U t'i+i,  we  can  rewrite  (2.87)  as 


t,+l 


S^'V1  - 2SJ  ' ' + |Vjd'z|  + |ViD't1+1| 


(2.88) 


Substituting  (2.86)  into  (2.88)  we  get 


5j  1 


t.+l 


4Sj"t1'1  + 3 | V j d ' Z I + 2S  j_t]_  + |ViD't1+1| 


(2.89) 
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4sj"tl-l  + 3 | V j d ' z I + Si^l+l. 


(2.90) 


In  subproblem  q^-1,  we  have 

+ (n-l)  |vjc'zl  + S^1-1 

where  n=  2<3q"tl-  Note  here  that  < 0. 

Substituting  (2.85)  into  (2.91),  we  get 

Sj"^-1  ^ - |Vjd'ZI  + s^1-1. 


(2.91) 


(2.92) 


In  subproblem  q^,  since  iG  subtree  Tj^  is  pruned  from  subtree 

I t I — 1 

Tj  91,  thus 


Sj  ql  ^ -|VjD  Zl . 


(2.93) 


Now,  starting  from  subproblem  qi,  and  using  the  same  argument  as  above, 
we  can  show  that  for  subproblems  q,  where  iGZj9,  we  have 


t t t f . , III  — , 

Sj  9 < - | VjD  z|. 

For  subproblems  t,  where  i^Zj1-,  we  have 


(2.94) 


f t i *-  , III  — , l 

Sj  C * -I VjD  z|  + Si 


:t 


(2.95) 


where  S^  < 0. 

Thus,  for  case  a),  where  iGZj T , we  can  use  (2.94)  by  letting 
r=q , and  have 

Sj  r ^ -|vjD  z|. 


(2.96) 
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Similarly,  for  case  b) , where  i^Zj r , we  can  use  (2.95)  by  letting 
t— t,  and  have 

Sj  r < - | Vj D z|  + Sir.  (2.97) 

Both  (2.96)  and  (2.97)  is  a contradiction  to  the  fact  that  all  Sj ' ' 
values  are  positive  in  the  original  problem.  This  completes  the  proof 
for  cases  1 and  2. 


Case  3 : Consider  a node  i€Tj  such  that  i^Zj2  but  ieZj T . (For  example 

node  3 in  Figure  2.15c). 

Case  4 : Consider  a node  ieTj ' such  that  i^Zjz  and  i^Zj r . (For  example 

node  4 in  Figure  2 . 15d) 

We  will  use  the  same  proof  for  cases  3 and  4. 

On  the  contrary,  we  will  assume  that,  Sj,z<  - | V J ' z | ' 

At  subproblem  z,  we  can  divide  Vjp’2  into  disjoint  sets  as 


vjo'Z  - vi2  U (vjD'Z\V-z). 


We  will  consider  the  following  subproblems,  qp  , tj. , q2  , t2  , . . . q^,  t^ 
where  z < qp  < t\  < . . . < t^  < r . In  subproblems  z to  qq-1,  we  have 
i^Zj  . In  subproblems  qp  to  tj_-l,  we  have  ieZj  . In  general,  in 
subproblems  q^  to  t^-l,  we  have  ieZj , and  in  subproblems  t^  to  q^+p  -1, 
we  have  i^Zj . 

Since  the  subtree  Tj  does  not  change  till  subproblem  q^ , we  can 
show  that  at  subproblem  qq-1,  we  have 


sj  ' S -|vji'2\vi*|  - |vi=|  . 


(2.98) 
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At  subproblem  , ieZjq^,  thus  > 0,  and  subtree  is  pruned  by 

procedure  PRUNE2 . Thus,  the  most  increase  in  Sj  is  achieved  when 

sj  " q1  = 2Sj  " cl1'1  + |Vjd'z\v[z|  + |v[z|  - s[ql 


^ -|VjD'zWiz|  - |viz| 


The  most  increase  in  Sj  is  achieved  when 


sj"ql+1  - 2s-"qi  + l vj d * i1+1 1 


2S-"ql  + |Vjd'z\v[z| 


* -|v-d'z\v[z|  - |viz| 


In  subproblem  t^-1,  we  have 


+ | Vj o* z\vj_! 


* -|vji'zwi*Mviz|. 


In  subproblem  t , subtree  T^t^  is  included  in  Tj  c . Thus  we  know  that 
Sj^  < 0,  and  the  most  increase  in  (Sj  c-*-)  is  achieved  by 


tl 


+ |Vji’z\vi*|  + si 


;ti 


* - |vjo'z\vizl  - |v[z|. 


Continuing  in  this  manner,  we  can  see  that  for  all  the  remaining 
problems  Sj  will  be  less  than  -|Vjd  z \ V ^ z | - |V^Z|  contradicting  to 

t t t _ 

the  fact  that  all  Sj  T values  are  positive  in  the  original  problem. 
This  completes  the  proof  of  the  lemma. 


69 


a)  ieZj z and  ieZj r 


b)  ieZjz  and  i^Zjr 


c)  i^ZjZ  and  ieZj T 


d)  i^Zj z and  i^Zj  r 


Figure  2.15  Subtrees  T j ' ' ' 


at  Subproblems  z and  r 
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Lemma  2 . 14  : At  the  beginning  of  subproblem  z,  (Smax,az)B  is 
bounded  by  2|VaD' | + |VaB' | . 

Proof  : Note  that  from  Equations  (2.79)  and  (2.80),  we  have 

(Smax,ocZ)B  - (Sa'z)B  - Z Sj"z  + | Zaz  | 

jeZaz 

< (Sa'z)B  + S |V-d'z|  + I Z<xz|  (from  Lemma  2.13) 

jeZ<xz 

- IvocdI  + 2 I Vj d z|  + | Zaz | (from  Corollary  2.1) 

jeZaz 

- I ^aD I + 2 |Vjd  z|  + |VaB|  (from  the  definition 

jeZCTz  of  the  set  Zaz) 

2!  2 Kd|  + |V«B|  .■ 

Lemma  2 . 15  Smaxz , total  deficiency  at  the  beginning  of  subproblem 

z,  z > 1 is  bounded  by  3 | V | . 

Proof:  Note  that  S^axz  - S (Smax>a)B 

cceQz 

For  i,jeQz,  we  have  n Vj  = 0.  Thus,  using  Lemma  2.14,  we  have 


BmaxZ  — 3 | V | . 
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Theorem  2.5,  Total  number  of  pivots  on  the  new  scaling  algorithm, 
DUALSC  is  bounded  by  3(r+l)|V|^.  Moreover,  the  algorithm  requires 
0[(r+l)|Vp]  arithmetic  operations. 

Proof:  Same  as  proof  of  Theorem  2.4. 

The  bound  we  have  obtained  for  the  algorithm  DUALSC  depends  on 
the  logarithm  of  the  total  supplies  (or  demands) , thus  DUALSC  is  not  a 
genuinely  polynomial  algorithm.  However,  when  we  compare  this  bound 
with  the  computational  bounds  of  the  other  dual  simplex  algorithms  in 
the  literature,  we  see  that  it  is  an  0(|V|^)  better  than  Ikura  and 
Nemhauser  s algorithm  (designed  for  transportation  problems  only) . 

Orlin  s two  genuinely  polynomial  dual  simplex  algorithms  have 
computational  bounds  of  0( |V^| log|V| ) and  0( | | log | V | ) when 
BIT(b)  —oo.  We  want  to  restate  here  that  BIT(b)  is  the  minimum  value  of 
s such  that  2sb  is  integer  valued,  and  b is  the  vector  of  supply  and 
demand  values  with  |bj_|  < 1.  We  believe  that  in  most  cases  BIT(b)  will 
be  equal  to  infinity.  Thus  algorithm  DUALSC  compares  very  favorably 
with  all  the  other  dual  simplex  algorithms.  To  compare  our  algorithm 
empirically,  we  conducted  a computational  experimentation  which 
compares  all  the  dual  simplex  algorithms  on  randomly  generated 
problems.  We  present  this  experimentation  in  Section  2.8. 
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2,7  A Heuristic  Leaving-Arc  Selection  Rule 

Before  conducting  our  computational  experience , we  introduce 
another  dual  simplex  algorithm  that  differs  from  DUAL  described  in 
Section  2.4.5  in  the  determination  of  the  leaving  arc. 

In  this  algorithm,  step  la.  of  algorithm  DUAL  is  changed  as 
follows  : 

Step  la  ) 

Let  Q = { i : i<=VF,  S±  > 0 } 

If  p'  =0  THEN 

Find  Tp  such  that  D(p)=max  { D ( j ) } . 

jeQ 

ELSE 

Find  Tp  such  that  D(p)=  max  { D ( j ) } 

jeQnPp' 

END  IF 

We  will  call  this  algorithm  MAXDEP  since  the  node  to  be  pivoted 
on  is  the  one  with  the  maximum  depth  among  all  the  possible  candidates 
This  selection  rule  corresponds  to  deleting  the  first  found  arc  with  a 
negative  flow.  Even  though  we  could  not  obtain  a polynomial  bound  for 
this  algorithm,  we  will  include  it  in  the  computational  comparisons 
since  MAXDEP  avoids  the  computation  of  the  Sj/  values  and  may  require 
less  computation  time. 
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2 . 8 . Computational  Experience 

Our  purpose  is  to  compare  the  algorithms  DUAL  and  DUALSC  given  in 
Sections  2.4.5  and  2.6  respectively  with  other  dual  simplex  algorithms 
available  in  the  literature  and  with  a primal  simplex  algorithm  as  well 
as  the  algorithm  MAXDEP.  Two  performance  measurements  will  be 
considered;  number  of  simplex  (primal  or  dual)  pivots  and  the  CPU  time. 

The  only  two  dual  simplex  algorithms  for  minimum  cost  flow  problem 
with  polynomial  bounds  in  the  literature  are  given  by  Orlin  [29],  We 
included  them  in  our  computational  comparison  and  called  them  0RLIN1 
and  0RLIN2.  The  primal  simplex  algorithm  that  we  included  in  our  study 
is  the  GNET  program  [16],  a readily  available  state-of-the-art  code. 

In  the  algorithms  DUAL  and  DUALSC,  we  used  three  node -length 
arrays  to  represent  the  spanning  tree.  These  are  the  predecessor,  the 
thread  and  the  last  node  in  the  subtree  indices.  Three  node -length 
arrays  are  used  to  store  the  dual  values,  and  the  S^  and  S^  values.  For 
the  algorithm  DUALSC,  two  more  node -length  arrays  are  used  to  store  the 
S^  and  Sj^  values  for  the  final  subproblem.  Three  arc -length  arrays  are 
used  to  store  the  input  data.  They  are  from-node,  to-node  and  cost.  Two 
node -length  arrays  are  used  as  pointers  to  from-node  and  to-node 
arrays.  Weights  of  nodes  are  stored  in  one  node-length  array.  We  also 
used  other  node-length  arrays  (to  store  the  inverse  thread,  a flag  to 
indicate  the  nodes  in  a given  subtree,  and  to  store  the  orientation  of 
the  arcs)  to  reduce  the  computational  time  further.  We  found  that  these 
additional  arrays  reduced  the  overall  computation  time  even  though 
extra  time  was  spent  to  update  them. 
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The  coding  of  0RLIN1  and  0RLIN2  are  straightforward 
implementations  of  the  algorithms.  In  these  algorithms,  the  spanning 
tree  is  represented  with  three  node -length  arrays;  the  predecessor,  the 
thread  and  the  depth  indices. 

We  now  summarize  the  memory  requirements  of  all  each  code. 

Table  2.1:  Memory  Requirements  for  each  of  the  Codes. 


Code 

Memory  Locations  Used 

DUAL 

13  |V| 

+ 3 |A  | 

DUALS C 

15  |V| 

+ 3|a| 

GNET 

9 |V| 

+ 3|A| 

MAXDEP 

12  |V| 

+ 3|A| 

0RLIN1 

13  |V| 

+ 3|A| 

0RLIN2 

13  |V| 

+ 3 |A  | 

All  of  the  algorithms  tested  are  coded  in  FORTRAN  77  and  are  run 
on  a VAX  11/750  using  the  VAX- 11  Fortran  V3.0  compiler.  All  runs  are 
made  overnight  when  we  were  the  only  users  of  the  system.  All 
algorithms  are  compared  on  the  same  set  of  test  problems.  Test  problems 
are  randomly  generated  by  the  NETGEN  program  developed  in  [36].  For 
minimum  cost  flow  problems,  the  input  to  NETGEN  problem  are  the 
following;  number  of  source  nodes,  number  of  demand  nodes,  number  of 
transshipment  source  nodes,  number  of  transshipment  demand  nodes, 
number  of  arcs,  cost  range  and  total  supply  (or  demand).  Here 
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transshipment  source  (demand)  nodes  are  those  source  (demand)  nodes 
into  (from)  which  arcs  are  allowed  to  lead  (emanate). 

The  computational  studies  are  carried  out  in  two  phases . In  the 
first  phase,  the  algorithms  are  compared  for  small  problems.  The 
characteristics  of  the  problems  generated  in  phase  one  are  given  in 
Table  2.2. 


Table  2,2:  Characteristics  of  the  Problems  Generated  in  Phase  1. 


Problem 

Type 

No . of 
Nodes 

No . of 
Source 
Nodes 

No . of 
Demand 
Nodes 

No . of 
Trans . 
Source 
Nodes 

No . of 
Trans . 
Demand 
Nodes 

No . of 
Arcs 

Total 

Supply 

Cost 

Range 

M-l 

50 

5 

5 

5 

5 

2000 

10000 

1-100 

M-2 

50 

5 

5 

5 

5 

1500 

10000 

1-100 

M-3 

50 

5 

5 

5 

5 

750 

10000 

1-100 

M-4 

75 

5 

9 

5 

9 

4500 

20000 

1-100 

M-5 

75 

5 

9 

5 

9 

3375 

20000 

1-100 

M-6 

75 

5 

9 

5 

9 

1688 

20000 

1-100 

M-7 

100 

10 

20 

5 

5 

8000 

30000 

1-100 

M-8 

100 

10 

20 

5 

5 

6000 

30000 

1-100 

M-9 

100 

10 

20 

5 

5 

3000 

30000 

1-100 

Note  that  problems  M-l,  M-4,  and  M-7  are  dense  problems  whereas 
M-3,  M-6  and  M-9  are  sparse  problems. 

Results  of  the  experiments  in  phase  one  are  summarized  in  Table 
2.3  where  the  mean  values  are  obtained  by  averaging  ten  runs. 

The  CPU  times  that  are  reported  are  in  seconds , and  they  do  not 
include  the  input-output  time.  The  time  statistics  are  collected  using 
the  library  routines  LIB$INIT_TIMER  and  LIB$STAT_TIMER . 
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Table  2.3.  The  Performance  of  Algorithms  in  Phase  1 


No. 

of  Pivots 

CPU  Times 

Min 

Max 

Mean 

Min 

Max 

Mean 

Problem  M-l 

DUAL 

12 

33 

23.8 

0.29 

0.48 

0.43 

DUALS C 

15 

34 

25.5 

0.35 

0.59 

0.52 

GNET 

63 

83 

74.3 

0.38 

0.58 

0.51 

MAXDEP 

20 

60 

42.3 

0.28 

0.51 

0.47 

0RLIN1 

14 

38 

25.7 

0.99 

1.28 

1.18 

0RLIN2 

16 

39 

26.1 

1.01 

1.32 

1.20 

Problem  M-2 

DUAL 

10 

30 

21.9 

0.20 

0.39 

0.31 

DUALS C 

10 

34 

22.4 

0.27 

0.57 

0.41 

GNET 

57 

91 

76.0 

0.28 

0.48 

0.40 

MAXDEP 

19 

42 

27.9 

0.21 

0.46 

0.33 

0RLIN1 

13 

35 

23.4 

0.99 

1.20 

1.11 

0RLIN2 

17 

36 

24.1 

1.01 

1.23 

1.14 

Problem  M-3 

DUAL 

10 

31 

18.4 

0.14 

0.27 

0.18 

DUALS C 

10 

33 

23.0 

0.21 

0.40 

0.31 

GNET 

65 

94 

75.4 

0.21 

0.31 

0.25 

MAXDEP 

11 

31 

23.8 

0.12 

0.25 

0.19 

0RLIN1 

10 

33 

19.2 

0.59 

0.73 

0.65 

0RLIN2 

13 

35 

21.2 

0.60 

0.77 

0.67 

Problem  M-4 

DUAL 

33 

43 

36.8 

0.87 

1.42 

1.11 

DUALS C 

32 

43 

39.2 

0.91 

1.62 

1.32 

GNET 

115 

146 

129.3 

1.15 

1.65 

1.39 

MAXDEP 

34 

56 

45.5 

0.79 

1.39 

1.09 

0RLIN1 

35 

45 

38.4 

2.98 

3.45 

3.12 

0RLIN2 

40 

48 

42.5 

3.05 

3.62 

3.27 
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Table  2.3.  (continued) 


No.  of  Pivots 
Min  Max  Mean 

CPU  Times 
Min  Max  Mean 

Problem  M-5 

DUAL 

25 

51 

35.0 

0.63 

1.03 

0.84 

DUALS C 

30 

42 

34.0 

0.69 

1.06 

0.92 

GNET 

112 

137 

123.6 

0.90 

1.15 

1.03 

MAXDEP 

36 

68 

47.6 

0.55 

1.18 

0.89 

0RLIN1 

27 

54 

38.7 

2.51 

3.03 

2.87 

ORLIN2 

30 

61 

42.1 

2.48 

3.13 

2.91 

Problem  M-6 

DUAL 

18 

34 

26.8 

0.38 

0.46 

0.43 

DUALS C 

21 

40 

33.5 

0.45 

0.87 

0.71 

GNET 

91 

131 

117.3 

0.45 

0.68 

0.58 

MAXDEP 

28 

34 

30.5 

0.37 

0.50 

0.44 

0RLIN1 

21 

39 

30.2 

1.44 

1.98 

1.70 

ORLIN2 

24 

45 

34.1 

1.65 

1.89 

1.79 

Problem  M-7 

DUAL 

40 

59 

49.8 

1.71 

2.69 

2.13 

DUALS C 

35 

60 

50.0 

1.94 

3.14 

2.66 

GNET 

150 

182 

168.5 

2.18 

2.99 

2.68 

MAXDEP 

41 

76 

61.5 

1.63 

3.33 

2.45 

0RLIN1 

42 

65 

52.1 

3.99 

5.03 

4.75 

ORLIN2 

46 

71 

58.0 

4.02 

5.14 

4.84 

Problem  M-8 

DUAL 

32 

69 

47.3 

1.30 

2.12 

1.60 

DUALS C 

39 

83 

59.1 

1.48 

3.47 

2.29 

GNET 

151 

176 

158.1 

1.74 

2.39 

2.01 

MAXDEP 

49 

130 

80.6 

1.25 

2.98 

2.17 

ORLIN1 

28 

74 

50.2 

2.94 

4.76 

3.91 

0RLIN2 

31 

77 

53.5 

2.97 

4.87 

3.93 

Problem  M-9 

DUAL 

28 

64 

44.4 

0.55 

1.28 

0.89 

DUALS C 

29 

68 

49.0 

0.77 

1.52 

1.13 

GNET 

147 

174 

164.6 

1.16 

1.27 

1.20 

MAXDEP 

43 

117 

68.4 

0.55 

1.62 

0.99 

ORLIN1 

30 

69 

49.2 

3.12 

3.99 

3.69 

ORLIN2 

32 

73 

51.2 

3.34 

3.98 

3.74 
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In  almost  all  but  one  of  the  nine  problem  types,  algorithm  DUAL 
was  faster  than  all  the  other  algorithms  tested.  When  we  examine  the 
ratios  of  CPU  times  of  other  algorithms  to  the  CPU  times  of  algorithm 
DUAL  in  Table  2.4,  we  see  that  DUAL  is  two  to  four  times  faster  than 
algorithms  0RLIN1  and  0RLIN2 . Algorithm  DUALSC  is  faster  in  dense 
problems  and  slower  in  sparse  problems.  Algorithm  MAXDEP  is  almost  as 
fast  as  algorithm  DUAL.  In  problem  type  M-4,  it  is  even  faster. 
However,  in  problem  type  M-8,  it  is  relatively  slower  than  DUAL.  In 
this  problem  type,  MAXDEP  required  almost  twice  the  pivots  algorithm 
DUAL  required. 

Next,  we  will  examine  the  CPU  time  per  pivot  for  each  of  the 
algorithms  in  Table  2.5. 

Algorithm  GNET  requires  three  to  four  times  more  pivots  compared 
to  the  other  algorithms.  However,  it  requires  the  minimum  amount  of 
time  per  pivot. 

Table  2.4,  Ratio  of  the  CPU  Times  of  Algorithms  to  the  CPU  Time  of 
Algorithm  DUAL 


Problem  Type 

DUALSC 

GNET 

MAXDEP 

0RLIN1 

0RLIN2 

M-l 

1.21 

1.19 

1.09 

2.74 

2.79 

M-2 

1.32 

1.29 

1.06 

3.58 

3.68 

M-3 

1.72 

1.38 

1.06 

3.61 

3.72 

M-4 

1.19 

1.25 

0.98 

2.81 

2.95 

M-5 

1.09 

1.23 

1.06 

3.41 

3.46 

M-6 

1.65 

1.35 

1.02 

3.95 

4.16 

M-7 

1.24 

1.26 

1.15 

2.22 

2.27 

M-8 

1.43 

1.26 

1.36 

3.21 

3.24 

M-9 

1.27 

1.35 

1.11 

4.15 

4.20 
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Table  2.5.  Time  per  Pivot  for  each  of  the  Algorithms 


Problem  Type 

DUAL 

DUALS C 

GNET 

MAXDEP 

0RLIN1 

0RLIN2 

M-l 

.018 

.021 

.007 

.011 

.046 

.046 

M-2 

.014 

.018 

.005 

.012 

.047 

.047 

M-3 

.010 

.013 

.003 

.008 

.033 

.031 

M-4 

.030 

.034 

.011 

.024 

.081 

.076 

M-5 

.024 

.027 

.008 

.019 

.074 

.069 

M-6 

.016 

.021 

.005 

.014 

.056 

.052 

M-7 

.043 

.053 

.016 

.039 

.085 

.083 

M-8 

.034 

.039 

.013 

.027 

.078 

.073 

M-9 

.020 

.023 

.007 

.014 

.075 

.073 

In  phase  two  of  the  experimentation,  we  will  compare  the 
algorithms  for  larger  problems.  However,  since  the  performance  of 
algorithms  0RLIN1  and  0RLIN2  are  very  similar  and  poorer  compared  to 
the  other  algorithms  for  small  problems,  we  will  include  only  the 
0RLIN1  algorithm  in  the  second  phase. 

In  phase  two,  we  generated  two  groups  of  test  problems,  one  dense 
and  the  other  sparse.  The  characteristics  of  these  problems  are 
summarized  in  Table  2.6. 


Table  2 . 6 . Characteristics  of  the  Problems  Generated  in  Phase  2 


Problem 

Type 

No . of 
Nodes 

No . of 
Source 
Nodes 

No . of 
Demand 
Nodes 

No . of 
Trans . 
Source 
Nodes 

No . of 
Trans . 
Demand 
Nodes 

No . of 
Arcs 

Total 

Supply 

Cost 

Range 

MCD-1 

150 

20 

40 

10 

20 

18000 

10000 

1-1000 

MCD-2 

150 

20 

50 

10 

30 

16000 

20000 

1-1000 

MCD-3 

200 

10 

50 

10 

30 

28000 

30000 

1-1000 

MCD-4 

200 

10 

20 

5 

20 

32000 

30000 

1-1000 

MCS-1 

1000 

70 

50 

30 

20 

10000 

10000 

1-1000 

MCS-2 

1500 

100 

100 

50 

50 

22500 

20000 

1-1000 

MCS-3 

2000 

200 

200 

90 

90 

30000 

30000 

1-1000 
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Problems  MCD-1  to  MCD-4  are  dense  problems  whereas  MCS-1  to  MCS-3  are 
sparse  problems.  Results  of  the  experiments  in  phase  two  are  summarized 
in  Table  2.7. 

Table  2.7.  Performance  of  Algorithms  in  Phase  Two 


Min 

No . of 
Max 

Pivots 

Mean 

Min 

CPU  Times 
Max 

Mean 

Problem  MCD-1 

DUAL 

89 

110 

96.7 

6.44 

9.15 

7.75 

DUALS C 

88 

101 

94.2 

6.28 

10.26 

7.84 

GNET 

244 

277 

259.0 

6.48 

7.28 

7.02 

MAXDEP 

128 

148 

139.5 

6.69 

8.35 

7.36 

0RLIN1 

95 

130 

102.5 

15.11 

24.71 

20.16 

Problem  MCD-2 

DUAL 

87 

115 

98.1 

6.56 

8.43 

7.12 

DUALS C 

86 

121 

103.4 

6.68 

9.38 

7.25 

GNET 

260 

294 

275.2 

6.15 

8.13 

6.99 

MAXDEP 

161 

185 

172.4 

6.99 

12.66 

9.65 

0RLIN1 

93 

116 

103.5 

13.42 

22.13 

19.15 

Problem  MCD-3 

DUAL 

46 

89 

72.1 

7.36 

13.15 

10.54 

DUALS C 

63 

110 

87.1 

9.09 

17.74 

13.23 

GNET 

306 

380 

344.5 

9.81 

12.59 

11.21 

MAXDEP 

90 

114 

107.3 

9.54 

15.88 

12.58 

0RLIN1 

51 

94 

73.2 

25.13 

32.71 

29.62 

Problem  MCD-4 

DUAL 

81 

143 

125.2 

12.15 

23.27 

16.11 

DUALS C 

98 

186 

138.1 

13.76 

32.39 

19.49 

GNET 

316 

390 

352.5 

11.47 

16.87 

14.12 

MAXDEP 

161 

292 

250.4 

14.67 

27.88 

21.23 

ORLIN1 

98 

146 

126.2 

30.12 

71.34 

53.24 
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Table  2.7  (continued) 


Min 

No . of 
Max 

Pivots 

Mean 

CPU 

Min 

Times 

Max 

Mean 

Problem 

MCS-1 

DUAL 

768 

968 

840.2 

44.07 

60.31 

54.12 

DUALS C 

804 

1248 

934.8 

63.76 

90.12 

75.34 

GNET 

2100 

2314 

2218.9 

24.16 

26.14 

25.76 

MAXDEP 

2661 

3212 

2812.5 

42.12 

60.12 

50.70 

0RLIN1 

2010 

2157 

2098.1 

224.60 

255.16 

250.91 

Problem 

MCS-2 

DUAL 

3100 

3321 

3214.9 

90.60 

162.44 

134.23 

DUALS C 

3134 

3414 

3297.1 

130.36 

175.45 

157.39 

GNET 

3301 

3547 

3421.6 

40.31 

46.12 

43.60 

MAXDEP 

4276 

5412 

4892.2 

75.45 

151.31 

110.25 

0RLIN1 

3623 

3817 

3756.1 

403.18 

501.90 

489.10 

Problem 

MCS-3 

DUAL 

4215 

4515 

4412.1 

200.45 

250.20 

244.17 

DUALS C 

4314 

4575 

4476.8 

278.56 

302.71 

295.16 

GNET 

4612 

4819 

4798.0 

80.16 

85.15 

83.16 

MAXDEP 

5513 

5812 

5676.2 

215.23 

254.65 

223.78 

0RLIN1 

4278 

4565 

4423.4 

696.52 

787.82 

723.02 

Within  the  size  range  of  problems  tested,  we  can  see  that  for  dense 
problems,  algorithms  DUAL,DUALSC,  MAXDEP  and  GNET  are  very  competitive. 

In  particular,  when  we  compare  algorithms  DUAL  and  GNET,  we  see  that  for 
problems  MCD-1,  MCD-2,  and  MCD-4,  GNET  is  slightly  faster  than  DUAL 
whereas  for  problem  type  MCD-3,  DUAL  is  slightly  faster  than  GNET.  When 
we  compare  them  in  terms  of  number  of  pivots  criterion,  we  see  that  GNET 
required  three  to  five  times  more  pivots  than  DUAL  for  dense 
problems .Algorithm  MAXDEP  performed  as  good  as  DUAL  and  GNET  in  the  first 
three  dense  problems. 
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Algorithm  0RLIN1  performed  poorer  compared  to  the  other  algorithms 
for  all  the  problem  types  in  phase  two.  It  required  three  to  five  times 
more  CPU  time. 

When  we  examine  the  performance  of  algorithms  for  sparse  problems, 
we  see  that  GNET  was  superior  compared  to  all  the  other  algorithms.  It  is 
twice  or  three  times  as  fast  as  algorithm  DUAL.  For  problem  types  MCS-2 
and  MCS-3,  the  number  of  pivots  GNET  required  and  the  number  of  pivots 
DUAL  required  are  close  to  each  other.  In  these  problems,  GNET  is  three 
times  faster  than  DUAL. 

Based  on  the  experimentation,  we  can  conclude  that  algorithm  DUAL  is 
faster  than  GNET  for  small  problems.  For  large,  dense  problems,  they  are 
competitive.  However,  for  large,  sparse  problems,  algorithm  GNET  is  two 
to  three  times  faster  than  algorithm  DUAL. 

Now,  we  will  compare  the  results  obtained  from  this  study  to  the 
results  given  by  Glover,  Karney  and  Klingman  [10].  In  their  study,  they 
compared  the  primal-simplex  algorithm  PNET  with  other  algorithms,  mainly 
with  a dual-simplex  algorithm  DNET,  and  out-of-kilter  algorithms.  Their 
conclusion  was  PNET  was  superior  to  all  the  other  algorithms.  In 
particular,  their  results  showed  that  PNET  is  6-12  times  faster  than 
DNET.  There  is  a major  difference  between  their  conclusion  and  ours  on 
performance  of  the  dual  codes.  The  basic  reasons  for  this  gap  can  be 
summarized  in  the  following  way  : 

a)  In  DNET,  the  leaving  arc  is  determined  by  the  first  negative 
criterion.  In  DUAL  and  DUALSC  , the  leaving  arc  selection  rules  are 


different . 
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b)  They  only  tested  for  sparse  problems. 

c)  As  they  have  concluded  in  their  study,  the  representation  of  the  basis 
in  DNET  is  not  the  most  efficient  representation  for  the  dual  approach. 
Using  the  inverse -thread  index  reduces  the  CPU  time  a great  deal  in  DUAL 
and  DUALSC. 


f 
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CHAPTER  3 

THE  TRANSPORTATION  PROBLEM 


3.1  Formulation  of  the  Problem 


The  transportation  problem  is  a special  case  of  the  minimum  cost 
flow  problem.  Consider  the  network  representation  G = (V,E)  rooted  at 
node  0 as  discussed  in  section  2.3  with  V-  Vg  U Vq.  If  there  exists  a 
direct  arc  from  ieVg  to  jeVp,  then  (i,j)eE.  For  each  arc  (i,j)eE,  we 
associate  a nonnegative  integer  cjj  which  corresponds  to  a unit  cost 
over  the  arc  (i,j).  For  each  node  ieVg,  we  associate  a positive  integer 
a^  designating  the  supply  at  that  point  and  for  each  node  ieVp , we 
associate  a positive  integer  bj  corresponding  to  the  demand  at  that 
point . 


Then,  we  can  formulate  the  transportation  problem  (TP)  exactly 
the  same  as  the  (MCP)  formulation. 

(TP)  Minimize  2 cijxij  (3.1) 

(i.j)eE 

subject  to 


2 Xjj  - 2 xij 

(i.j)ex(i)  (i,j)ei9(i) 


-bH 


if  ieVs 
if  ieVD 


(3.2) 


(i,j)eE 


(3.3) 


xij  * 0 
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3,2  The  Dual  Simplex  Algorithms  for  the  Transportation  Problem 

The  steps  of  algorithms  DUAL,  described  in  Section  2.5.5  and 
DUALSC,  described  in  Section  2.6  can  be  applied  directly  for  the 
transportation  problem.  However,  the  construction  of  the  initial  tree 
can  be  changed  to  take  advantage  of  the  special  structure  of  the 
transportation  problem.  In  this  section,  we  will  only  restate  the 
construction  of  the  initial  tree,  that  is  the  initialization  step. 

t 

Let  Tq  be  a shortest  path  tree  of  G rooted  at  node  0 using  costs 
as  the  weights.  Such  a tree  can  easily  be  constructed  by  attaching  each 
supply  node  to  node  0 and  attaching  each  demand  node  j to  a supply  node 
x with  min{cjj  : ieVg)  - cxj • Let  w^  represent  the  negative  of  the 
shortest  path  from  node  0 to  node  i to  Tq  . Now  attach  each  source  node 
i,  which  is  not  attached  to  a demand  node,  to  a demand  node  k (after 
deleting  arc  ( 0 , i ) > with  w^  + c^  =■  min  { wj  + c^j  : jgVd)  and  set 
wi  “ wk  + cik- 

3.3  Computational  Experience 

In  this  section,  we  will  compare  algorithms  DUAL  and  DUALSC  with 
the  algorithms  discussed  in  section  2.7  on  transportation  problems. 
Moreover,  we  will  also  include  the  unsealed  and  scaled  versions,  named 
DTRAN  and  DTRANSC  respectively,  of  the  dual  simplex  algorithm  developed 
by  Ikura  and  Nemhauser  [32,33].  The  codings  of  DTRAN  and  DTRANSC  are 
the  original  ones  obtained  from  the  autors . 

Again,  all  of  the  algorithms  tested  are  coded  in  FORTRAN  77  and 
are  run  on  a VAX  11/750  using  the  VAX- 11  Fortran  V3.0  compiler.  All 
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runs  are  made  overnight  when  we  were  the  only  users  of  the  system.  All 
algorithms  are  compared  on  the  same  set  of  test  problems.  Test  problems 
are  generated  by  the  NETGEN  program. 

The  computational  studies  are  carried  out  in  two  phases.  In  the 
first  phase,  the  algorithms  are  compared  for  small  sized  problems.  The 
characteristics  of  the  problems  generated  in  phase  one  are  given  in 
Table  3.1. 

Table  3 . 1 Characteristics  of  the  Problems  Generated  in  Phase  One 


Problem 

Type 

Number  of 
Source 
Nodes 

Number  of 
Demand 
Nodes 

Number  of 
Arcs 

Total 

Supply 

Cost 

Range 

T-l 

70 

70 

4900 

10000 

1-100 

T-2 

70 

70 

2940 

10000 

1-100 

T-3 

70 

70 

1470 

10000 

1-100 

T-4 

50 

90 

4500 

10000 

1-100 

T-5 

50 

90 

2700 

10000 

1-100 

T-6 

50 

90 

1350 

10000 

1-100 

T-7 

90 

50 

4500 

10000 

1-100 

T-8 

90 

50 

2700 

10000 

1-100 

T-9 

90 

50 

1350 

10000 

1-100 

Note  that  problems  T-l,  T-4,  and  T-7  are  dense  problems  in  which 
the  density  equals  one.  Here  we  define  the  density  as  the  ratio  of  the 
number  of  actual  arcs  in  the  transportation  network  to  the  maximum 
possible  number  of  arcs.  Problems  T-3,  T-6  and  T-9  are  sparse  problems  in 
which  the  density  equals  0.3.  For  problems  T-2,  T-5  and  T-8,  the  density 
equals  0.6. 

Results  of  the  experiments  in  phase  one  are  summarized  in  Table  2.3 
where  the  mean  values  are  obtained  by  averaging  ten  runs. 

The  CPU  times  that  are  reported  are  in  seconds , and  they  do  not 
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include  the  input-output  time.  The  time  statistics  are  collected  using 
the  library  routines  LIB$INIT_TIMER  and  LIB$STAT_TIMER. 


Table  3.2  The  Performance  of  the  Algorithms  in  Phase  One. 


No. 

Min 

of  Pivots 
Max 

Mean 

CPU 

Min 

Times 

Max 

Mean 

Problem  T-l 

DUAL 

166 

223 

189.8 

2.03 

3.54 

2.84 

DUALS C 

172 

224 

199.8 

3.43 

5.42 

4.32 

DTRAN 

198 

250 

224.0 

3.54 

5.31 

4.42 

DTRANSC 

245 

276 

254.0 

3.72 

6.28 

4.58 

0RLIN1 

190 

223 

207.0 

6.01 

7.44 

6.55 

0RLIN2 

214 

259 

230.0 

6.09 

8.34 

7.33 

GNET 

231 

250 

239.0 

3.01 

3.46 

3.23 

MAXDEP 

166 

358 

251.6 

2.30 

3.94 

2.86 

Problem  T-2 

DUAL 

173 

189 

180.8 

1.52 

2.04 

1.77 

DUALS C 

186 

209 

196.8 

2.67 

3.77 

3.05 

DTRAN 

213 

257 

233.6 

2.91 

4.77 

3.72 

DTRANSC 

250 

280 

263.6 

2.84 

3.43 

3.20 

0RLIN1 

206 

245 

227.8 

6.23 

7.82 

6.73 

0RLIN2 

223 

286 

257.8 

6.53 

8.51 

7.21 

GNET 

251 

292 

275.2 

2.32 

2.86 

2.60 

MAXDEP 

199 

310 

247.4 

1.61 

2.61 

1.92 

Problem  T-3 

DUAL 

172 

190 

182.0 

1.08 

1.33 

1.21 

DUALS C 

177 

212 

189.8 

1.93 

2.35 

2.07 

DTRAN 

212 

237 

224.3 

1.93 

2.72 

2.23 

DTRANSC 

226 

280 

252.0 

1.52 

2.50 

2.01 

ORLIN1 

185 

247 

217.5 

4.19 

6.65 

5.22 

ORLIN2 

216 

282 

246.0 

4.49 

6.46 

5.44 

GNET 

239 

273 

250.5 

1.41 

1.63 

1.54 

MAXDEP 

231 

254 

239.3 

1.03 

1.27 

1.15 

Problem  T-4 

DUAL 

135 

173 

157.0 

1.79 

2.99 

2.48 

DUALS C 

139 

187 

158.0 

2.64 

3.97 

3.17 

DTRAN 

211 

229 

220.4 

3.32 

4.29 

3.87 

DTRANSC 

239 

270 

254.8 

3.41 

5.88 

4.17 

0RLIN1 

166 

200 

178.0 

5.16 

7.00 

6.02 

0RLIN2 

183 

236 

202.3 

5.84 

7.45 

6.38 

GNET 

234 

255 

241.4 

2.79 

3.23 

2.98 

MAXDEP 

197 

287 

213.8 

1.81 

2.69 

2.27 

Table  3.2  (continued) 


No. 

Min 

of  Pivots 
Max 

Mean 

CPU 

Min 

Times 

Max 

Mean 

Problem  T-5 

DUAL 

142 

191 

170.4 

1.27 

1.96 

1.76 

DUALS C 

171 

210 

182.3 

2.48 

2.90 

2.78 

DTRAN 

201 

257 

235.0 

2.53 

4.74 

3.87 

DTRANSC 

206 

255 

228.4 

1.67 

2.98 

2.57 

ORLIN1 

209 

262 

234.6 

6.49 

7.88 

6.94 

ORLIN2 

228 

298 

258.2 

5.78 

8.03 

7.04 

GNET 

219 

295 

260.8 

2.14 

2.62 

2.40 

MAXDEP 

173 

287 

234.0 

1.06 

2.24 

1.76 

Problem  T-6 

DUAL 

141 

183 

166.6 

0.98 

1.35 

1.13 

DUALS C 

174 

191 

181.6 

1.87 

2.21 

1.62 

DTRAN 

235 

260 

246.0 

2.23 

3.16 

2.84 

DTRANSC 

211 

270 

248.8 

1.59 

2.69 

2.18 

ORLIN1 

189 

211 

202.8 

4.53 

5.21 

4.86 

ORLIN2 

208 

255 

235.3 

4.45 

5.66 

4.99 

GNET 

249 

281 

270.0 

1.52 

1.75 

1.63 

MAXDEP 

180 

283 

232.0 

0.83 

1.45 

1.10 

Problem  T-7 

DUAL 

174 

204 

192.8 

2.04 

2.90 

2.50 

DUALS C 

181 

216 

197.4 

2.75 

4.02 

3.49 

DTRAN 

204 

234 

221.0 

3.27 

5.11 

4.05 

DTRANSC 

222 

246 

231.6 

2.43 

4.14 

3.52 

ORLIN1 

183 

214 

203.2 

5.39 

6.87 

6.28 

ORLIN2 

187 

212 

204.5 

5.41 

6.72 

6.23 

GNET 

224 

286 

244.8 

2.87 

3.99 

3.35 

MAXDEP 

240 

288 

262.2 

1.85 

3.35 

2.77 

Problem  T-8 

DUAL 

163 

219 

187.0 

1.18 

2.69 

1.91 

DUALS C 

165 

224 

191.5 

2.40 

3.08 

2.66 

DTRAN 

198 

229 

213.0 

2.65 

3.41 

2.96 

DTRANSC 

211 

232 

221.8 

1.80 

3.16 

2.52 

ORLIN1 

197 

232 

226.0 

5.97 

8.23 

6.84 

ORLIN2 

227 

283 

258.8 

6.21 

8.08 

6.88 

GNET 

242 

265 

250.5 

2.34 

2.56 

2.46 

MAXDEP 

195 

314 

254.5 

1.26 

3.04 

2.15 
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Table  3,2.  (continued) 


No. 

Min 

of  Pivots 
Max 

Mean 

CPU 

Min 

Times 

Max 

Mean 

Problem  T-9 

DUAL 

160 

199 

177.8 

0.84 

1.38 

1.14 

DUALSC 

166 

206 

181.3 

1.42 

2.13 

1.91 

DTRAN 

197 

208 

202.0 

1.61 

1.82 

1.71 

DTRANSC 

211 

234 

227.0 

1.44 

2.04 

1.75 

0RLIN1 

190 

258 

223.3 

4.39 

5.96 

5.30 

0RLIN2 

227 

280 

259.0 

4.09 

7.46 

5.97 

GNET 

216 

251 

233.0 

1.38 

1.64 

1.50 

MAXDEP 

204 

294 

241.8 

0.90 

1.79 

1.29 

In  all  but  three  of  the  nine  problem  types , algorithm  DUAL  was 

faster  than  all  the  other  algorithms  tested.  When  we  examine  the  ratios 

of  CPU  times  of  other  algorithms  to  the  CPU  times  of  algorithm  DUAL  in 

Table  3.3,  we  see  that  DUAL  is  three  to  five  times  faster  than  algorithms 

0RLIN1  and  0RLIN2 . Algorithm  DUALSC  is  faster  in  dense  problems  and 

slower  in  sparse  problems  whereas  GNET  is  faster  in  sparse  problems  and 

slower  in  denser  problems.  Algorithm  MAXDEP  is  as  fast  as  DUAL  and  in 

some  cases  faster  (Problems  M-5,  M-6,  and  M-9). 

When  we  compare  algorithms  DTRAN  and  DTRANSC,  we  see  that  DTRANSC 

runs  faster  than  DTRAN  in  most  of  the  problem  types  even  though  it 

required  more  iterations.  This  observance  agrees  with  the  conclusion  in 

[33].  When  we  compare  DUALSC  and  DTRANSC  we  see  that  there  is  no  major 

dominance  of  one  over  the  other.  In  some  problems  DUALSC  performed 

/ 

better,  and  in  others  DTRANSC  s performance  was  better.  However,  there  is 
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a observable  dominance  of  DUALSC  over  DTRAN.  DUALSC  is  20-50  percent 


faster  than  DTRAN. 


Table  3.3.  Ratio  of  the  CPU  Times  of  Algorithms  to  the  CPU  Time  of 
Algorithm  DUAL. 


DUALSC 

DTRAN 

DTRANSC 

0RLIN1 

0RLIN2 

GNET 

MAXDEP 

T-l 

1.52 

1.55 

1.61 

2.31 

2.58 

1.14 

1.01 

T-2 

1.72 

2.10 

1.81 

3.80 

4.07 

1.47 

1.08 

T-3 

1.71 

1.84 

1.66 

4.31 

4.50 

1.27 

0.95 

T-4 

1.28 

1.56 

1.68 

2.43 

2.57 

1.20 

0.91 

T-5 

1.57 

2.20 

1.46 

3.94 

4.00 

1.36 

1.00 

T-6 

1.43 

2.51 

1.93 

4.30 

4.42 

1.44 

0.97 

T-7 

1.40 

1.62 

1.41 

2.51 

2.75 

1.34 

1.11 

T-8 

1.39 

1.55 

1.32 

3.59 

3.60 

1.29 

1.13 

T-9 

1.57 

1.50 

1.53 

4.65 

5.24 

1.32 

1.13 

Next,  we  will  examine  the  CPU  time  per  pivot  for  each  of  the 
algorithms  in  Table  3.4. 

Algorithms  DUAL,  MAXDEP,  and  GNET  require  the  minimum  amount  of 

time  per  pivot.  Algorithm  DUALSC  requires  almost  twice  the  time  per  pivot 

of  algorithm  DUAL.  We  could  expect  this  since  in  algorithm  DUALSC,  for 

/ 

determining  the  leaving  arc,  we  need  to  determine  Sj^  values  for  the 
current  subproblem  and  the  final  subproblem. 


Table  3,4.  Time  per  Pivot  for  each  of  the  Algorithms 


DUAL 

DUALSC 

DTRAN 

DTRANSC 

0RLIN1 

0RLIN2 

GNET 

MAXDEP 

T-l 

0.015 

0.021 

0.019 

0.018 

0.031 

0.032 

0.013 

0.011 

T-2 

0.010 

0.015 

0.016 

0.012 

0.030 

0.028 

0.009 

0.008 

T-3 

0.007 

0.011 

0.010 

0.008 

0.024 

0.022 

0.006 

0.005 

T-4 

0.016 

0.020 

0.018 

0.016 

0.034 

0.032 

0.012 

0.011 

T-5 

0.010 

0.015 

0.016 

0.011 

0.030 

0.027 

0.009 

0.008 

T-6 

0.007 

0.009 

0.012 

0.009 

0.024 

0.021 

0.006 

0.005 

T-7 

0.013 

0.018 

0.018 

0.015 

0.033 

0.031 

0.014 

0.011 

T-8 

0.010 

0.014 

0.014 

0.011 

0.030 

0.027 

0.010 

0.008 

T-9 

0.006 

0.011 

0.008 

0.008 

0.024 

0.023 

0.006 

0.005 
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In  phase  two  of  the  experimentation,  we  will  compare  the  algorithms 
for  larger  problems.  However,  since  algorithms  0RLIN1  and  0RLIN2 
performed  poorer  compared  to  the  others  in  small  sized  problems,  we  will 
include  only  0RLIN1  in  the  second  phase.  In  phase  two,  we  generated  two 
groups  of  test  problems,  one  dense  and  the  other  sparse.  The 
characteristics  of  these  problems  are  summarizes  in  Table  3.5. 


Table  3 . 5 . Characteristics  of  the  Problems  Generated  in  Phase  2. 


Problem 

Type 

Number  of 
Source 
Nodes 

Number  of 
Demand 
Nodes 

Number  of 
Arcs 

Total 

Supply 

Cost  Range 

TD-1 

100 

200 

16000 

10000 

1-1000 

TD-2 

150 

150 

18000 

20000 

1-1000 

TD-3 

150 

200 

24000 

25000 

1-1000 

TD-4 

200 

200 

32000 

30000 

1-1000 

TS-1 

150 

150 

2250 

10000 

1-1000 

TS-2 

500 

1000 

5000 

10000 

1-1000 

TS-3 

1000 

1000 

10000 

20000 

1-1000 

TS-4 

1500 

1000 

15000 

30000 

1-1000 

Results  of  the  experiments  in  phase  two  are  summarized  in  Table  3.6. 


Table  3,6.  The  Performance  of  the  Algorithms  in  Phase  Two. 


No. 

Min 

of  Pivots 
Max 

Mean 

CPU  Times 
Min  Max 

Mean 

Problem 

TD-1 

DUAL 

369 

451 

407.0 

12.68 

17.87 

15.10 

DUALS C 

377 

449 

420.3 

19.24 

23.56 

20.71 

MAXDEP 

567 

743 

639.0 

14.83 

19.18 

16.76 

DTRANSC 

470 

646 

556.5 

21.89 

39.58 

30.71 

GNET 

549 

717 

618.5 

12.37 

16.50 

14.17 

0RLIN1 

381 

492 

455.1 

31:12 

42.65 

38.72 

Problem 

TD-2 

DUAL 

380 

532 

433.0 

13.16 

19.67 

15.63 

DUALS C 

387 

547 

469.1 

18.68 

25.21 

20.81 

MAXDEP 

570 

722 

652.0 

14.73 

22.06 

17.60 

DTRANSC 

507 

627 

554.4 

18.87 

32.14 

25.71 

GNET 

527 

630 

604.2 

13.54 

17.36 

15.84 

0RLIN1 

411 

562 

477.1 

42.16 

50.43 

47.90 

Problem 

TD-3 

DUAL 

462 

545 

514.8 

20.04 

27.38 

23.53 

DUALS C 

476 

548 

520.5 

39.12 

44.39 

40.29 

MAXDEP 

640 

796 

727.2 

19.49 

31.13 

24.38 

DTRANSC 

688 

851 

735.2 

49.78 

120.64 

74.25 

GNET 

712 

851 

761.5 

21.22 

27.02 

23.66 

ORLIN1 

500 

592 

567.0 

80.12 

91.78 

88.31 

Problem 

TD-4 

DUAL 

554 

631 

600.1 

28.71 

38.16 

34.91 

DUALS C 

579 

681 

652.7 

49.27 

64.12 

59.39 

MAXDEP 

799 

1797 

1193.4 

36.94 

132.40 

64.51 

DTRANSC 

744 

844 

804.6 

64.11 

119.72 

86.80 

GNET 

839 

928 

870.1 

31.92 

34.71 

32.12 

0RLIN1 

636 

712 

685.1 

91.34 

102.01 

98.91 
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Table  3 . 6 (continued) 


No. 

Min 

of  Pivots 
Max 

Mean 

CPU  Times 
Min  Max 

Mean 

Problem 

TS-1 

DUAL 

363 

477 

410.6 

3.25 

4.30 

3.63 

DUALS C 

401 

490 

450.2 

4.34 

5.28 

4.78 

MAXDEP 

506 

676 

570.6 

2.82 

4.02 

3.55 

DTRANSC 

535 

597 

568.1 

7.19 

7.87 

7.43 

GNET 

515 

624 

554.0 

3.29 

4.01 

3.59 

0RLIN1 

412 

499 

465.1 

8.11 

8.98 

8.65 

Problem 

TS-2 

DUAL 

1825 

1965 

1894.3 

43.22 

58.07 

50.60 

DUALS C 

1901 

2003 

1957.2 

68.23 

79.49 

70.32 

MAXDEP 

2661 

3030 

2406.9 

36.54 

55.69 

44.64 

DTRANSC 

2022 

2323 

2196.1 

65.10 

81.45 

74.80 

GNET 

2173 

2319 

2210.1 

22.08 

25.17 

23.43 

0RLIN1 

2145 

2218 

2193.0 

214.19 

245.90 

224.18 

Problem 

TS-3 

DUAL 

3126 

3328 

3180.5 

90.60 

122.92 

105.51 

DUALS C 

3176 

3341 

3190.1 

150.31 

200.35 

178.31 

MAXDEP 

4260 

5646 

4892.5 

75.88 

151.30 

106.33 

DTRANSC 

3233 

3612 

3451.3 

165.15 

220.27 

195.23 

GNET 

3305 

3649 

3449.1 

43.31 

47.40 

45.37 

0RLIN1 

3712 

3908 

3892.2 

411.55 

499.43 

455.89 

Problem 

TS-4 

DUAL 

4167 

4452 

4392.1 

170.12 

190.23 

183.26 

DUALS C 

4198 

4513 

4476.9 

240.54 

290.53 

275.39 

MAXDEP 

5412 

6587 

5912.7 

160.40 

195.76 

180.65 

DTRANSC 

4243 

4816 

4698.3 

280.32 

320.23 

302.88 

GNET 

4456 

4787 

4648.1 

70.25 

78.65 

72.88 

0RLIN1 

4423 

4654 

4532.3 

675.12 

712.48 

699.89 

When  we  examine  Table  3.6,  we  can  see  that  for  dense  problems, 


algorithms  DUAL  and  GNET  are  very  competitive  within  the  size  range  of 
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problems  tested.  Algorithm  MAXDEP  performed  as  good  as  DUAL  and  GNET  in 
the  first  three  dense  problems,  however,  for  problem  T-4,  DUAL  and  GNET 
were  almost  twice  as  fast  as  MAXDEP.  When  we  compare  algorithms  DUAL  and 
0RLIN1,  we  can  see  that  DUAL  is  almost  three  times  faster  than  0RLIN1  for 
dense  problems.  Also,  DUAL  is  almost  twice  as  fast  as  DTRANSC. 

When  we  examine  the  performance  of  algorithms  for  sparse  problems, 
we  see  that  GNET  was  superior  compared  to  all  the  other  algorithms  for 
larger  problems.  It  is  twice  as  fast  as  algorithm  DUAL.  However,  for 
problem  TS-1,  again  DUAL  and  GNET  are  competitive. 

In  conlusion,  algorithm  DUAL  is  faster  than  all  the  other  dual 
simplex  algorithms  in  the  literature  for  all  the  problem  types 
considered.  It  is  competitive  with  GNET  for  large  dense  problems, 
however,  for  large  sparse  problems  GNET  is  2-2.5  times  faster  than  DUAL. 

These  results  agree  with  our  conclusion  for  the  minimum  cost  flow 
problems  given  in  Section  2.8. 


CHAPTER  4 

THE  ASSIGNMENT  PROBLEM 


Assignment  problems  arise  when  a number  of  (jobs)  should  be 
assigned  to  a number  of  (machines)  with  the  minimum  cost  of  assignment. 

4.1.  The  Problem  Formulation 

The  assignment  problem  is  a special  case  of  the  transportation 
problem.  In  the  (TP)  formulation  in  Section  3.1,  if  - 1 for  all 
ieVg,  and  bj  - 1,  for  all  jeVj),  and  |Vg|  = |Vq|  = N,  then  the  problem 
is  called  an  assignment  problem  (AP) . 

4,2  Literature  Review 

For  assignment  problems,  several  algorithms  with  polynomial 
bounds  have  been  developed. 

First  algorithm  developed  for  assignment  problems  is  Kuhn  s [37] 
"Hungarian  Method"  with  a computational  bound  of  O(N^)  operations. 

Barr,  Glover  and  Klingman's  [20]  (AP-AB)  (Alternating  Path  - 
Alternating  Basis)  primal  simplex  algorithm  specialized  for  assignment 
problems  gave  good  results  computationally,  but  no  theoretical  bound  is 
known  for  this  algorithm. 
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Hung  and  Rom  [38]  developed  a recursive  algorithm  based  on 
relaxing  the  given  problem  which  exploits  the  (AP  - AB)  structure.  The 
algorithm  has  a computational  bound  of  0(N3)  operations. 

t t 

Engquist  s [39]  algorithm,  which  is  similar  to  Hung  and  Rom  s 
algorithm,  solves  the  assignment  problem  by  solving  successive  shortest 
path  problems.  The  bound  of  this  algorithm  is  0(N3)  operations. 

Hung  [40]  developed  a primal  simplex  algorithm  for  the  assignment 
problem  which  generates  at  most  N3ln  A bases  where  A is  the  difference 
in  the  objective  value  between  an  initial  extreme  point  and  the  optimal 
extreme  point.  While  preserving  (AP  - AB)  structure,  the  algorithm 
provides  a polynomial  bound  in  the  problem  size  and  the  encoding  of  the 
problem  data  by  applying  the  "Modified  Row  Most  Negative"  rule  in 
selecting  the  entering  arc. 

Bertsekas  s [41]  algorithm  for  the  assignment  problem  is  a 
primal-dual  algorithm  which  he  calls  as  an  "outpricing"  algorithm  with 
a complexity  of  0(N3)  + 0(RN3)  where  R is  an  integer  valued  number 
greater  than  all  cost  entries.  For  large  R,  the  algorithm  is  modified 
by  combining  with  the  Hungarian  method  to  give  a worst  case  complexity 
of  0(V3) . 

In  1985,  Balinski  [42]  provided  a dual  simplex  algorithm  for  the 
assignment  problem  that  terminates  in  at  most  (N-l)(N-2)/2  pivots.  The 
algorithm  is  guided  by  the  signatures,  where  signature  of  a tree  T is 
the  vector  of  its  source  node  degrees.  The  objective  of  the  method  is 
to  find  a tree  whose  signature  contains  exactly  one  1 and  otherwise 
2 s.  The  worst  case  bound  of  this  algorithm  agrees  with  the  bound  of 
our  algorithm.  The  main  difference  between  Balinski  s algorithm  and  the 


97 


DUAL  algorithm  is  that  in  his  algorithm,  one  can  select  a subtree  T-^, 
ieVs  for  pivoting  while  s[=0.  Algorithm  DUAL  requires  > 0 for  to 
be  selected  for  pivoting.  Due  to  this  difference,  Balinski's  algorithm 
cannot  detect  a dual  feasible  basis  which  admits  an  optimal  assignment. 
However,  in  1986,  Balinski  [43]  provided  a competitive  dual  simplex 
algorithm  with  the  same  computational  bound  as  the  one  in  [42] , which 
requires  to  be  greater  than  zero  if  is  to  be  selected  for 
pivoting.  Hence,  this  new  algorithm  is  equivalent  to  DUAL  applied  to 
assignment  problems.  However,  Balinski’s  algorithm  cannot  be  extended 
to  solve  the  minimum  cost  flow  problem  or  the  transportation  problem. 

4.3  A Dual  Simplex  Algorithm  for  the  Assignment  Problem  and  its 

Computational  Bound 

The  dual  simplex  algorithm  is  exactly  the  same  as  the  one  given 
in  Section  3.2. 

The  initial  total  deficiency  as  well  as  total  deficiency  during 
the  course  of  the  algorithm  is  defined  in  the  same  manner  as  in  the 
minimum  cost  flow  problem,  i.e. 

smax  " 2 (3.1) 

ieQ 

Note  that  Smax  is  equivalent  to  the  number  of  unsuccessful  assignments 
in  T , and  0 < Smax  < N - 2 . 

The  following  theorem  gives  the  upper  bound  on  the  number  of  dual 
simplex  pivots.  Let  u be  the  number  of  source  nodes  i in  T*  , the 
optimal  tree, such  that  (0,i)£  T*.  Let  Smax  be  the  initial  total 
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deficiency  as  defined  before,  and  q=|{i  : ieVp®,  Sj.  > 0)|,  that  is  q is 
the  number  of  nodes  i in  T®  with  Si>  0. 

Theorem  4.3.  The  total  number  of  dual  pivots  in  the  algorithm  for 
the  assignment  problem  is  bounded  by 

^max  + (u-1) (u-2)/2  - (q-l)(q-2)/2  (4.2) 

Proof:  At  the  beginning  of  the  algorithm,  there  can  be  at  most 
(q-1)  pivots  before  decreasing  Smax  by  one.  Thus  after  Smax  decreases 
by  one,  we  have  at  most  q+1  nodes  in  Vp  with  S^>0.  Thus  before  the  next 
reduction  occurs  in  Smax,  we  can  have  at  most  q pivots.  Continuing  with 
this  argument,  when  Smax  = 1 , we  have  at  least  two  source  nodes  in  Tp, 
thus  we  can  have  at  most  u-2  pivots  before  Smax  reduces  to  zero. 

Summing  up  all  the  pivots,  we  get 

u-2 

Smax  + 2 j - Smax  + (u-1) (u-2)/2  - (q-l)(q-2)/2  . ■ (4.3) 

j=q-i 

Lemma  4,3.  The  worst  case  bound  of  the  algorithm  is  (N-l)(N-2)/2 
pivots . 

Proof : Note  that  Smax  < N-2,  q > 1,  and  u < N-l.  Thus  by  letting 
Smax  -N-2,  q=l  and  u = N-l,  we  get 

(N-2)  + (N-3) (N-2)/2  = (N-l)(N-2)/2  . ■ 

4.4  Computational  Experience 

For  assignment  problems,  we  will  compare  the  algorithm  DUAL  with 
three  other  algorithms.  A primal  simplex  method,  GNET,  a primal-dual 
approach  coded  by  Burkard  and  Derigs  [44] , and  a dual  simplex 
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We  will  call  the  primal-dual  approach  algorithm  ASSIGN. 

Test  problems  are  randomly  generated  by  the  NETGEN  program 
developed  in  [36].  For  a fixed  number  of  nodes,  two  sets  of  problems 
are  generated;  one  being  a sparse  problem,  designated  by  AS  and  the 
other  being  denser,  designated  by  AD.  Characteristics  of  the  generated 
problems  are  given  in  Table  4.1. 


Table  4,1  Characteristics  of  the  Assignment  Problems  Generated 


Problem  Type 

No.  of  Nodes 

No.  of  Arcs 

Cost  Range 

AS  - 1 

100 

250 

1-100 

AS -2 

200 

1000 

1-100 

AS -3 

300 

2250 

1-100 

AS -4 

400 

4000 

1-100 

AD-1 

100 

1250 

1-100 

AD- 2 

200 

5000 

1-100 

AD- 3 

300 

11250 

1-100 

AD-4 

400 

20000 

1-100 

The  performance  of  the  algorithms  are  summarized  in  Table  4.2, 
where  the  mean  values  are  obtained  by  averaging  ten  runs.  Note  that  for 
algorithm  ASSIGN,  the  number  of  pivots  criterion  is  not  relevant  since 
it  does  not  perform  simplex- like  pivots. 

The  CPU  times  that  are  reported  are  in  seconds,  and  they  do  not 
include  the  input-output  time.  The  time  statistics  are  collected  using 
the  library  routines  LIB$INIT_TIMER  and  LI B $ S TAT_T IMER . 


Table  4.2  Performance  of  Algorithms  for  Assignment  Problems 


No 

Min 

of 

Max 

Pivots 

Mean 

Min 

CPU  Times 
Max 

Mean 

Problem  AS-1 

DUAL 

39 

73 

57.0 

0.10 

0.26 

0.18 

DTRAN 

92 

117 

104.6 

0.24 

0.43 

0.31 

ASSIGN 

* 

* 

* 

0.18 

0.26 

0.22 

GNET 

152 

188 

167.2 

0.37 

0.48 

0.43 

Problem  AS -2 

DUAL 

132 

180 

152.0 

0.73 

1.11 

0.92 

DTRAN 

222 

256 

240.2 

0.99 

1.91 

1.54 

SSIGN 

* 

* 

* 

0.99 

1.59 

1.29 

GNET 

335 

437 

384.4 

1.29 

1.79 

1.55 

Problem  AS -3 

DUAL 

215 

274 

237.6 

1.84 

2.49 

2.19 

DTRAN 

357 

402 

379.2 

2.97 

5.79 

4.71 

ASSIGN 

* 

* 

* 

2.94 

4.55 

3.24 

GNET 

523 

726 

637.6 

2.85 

3.95 

3.45 

Problem  AS -4 

DUAL 

296 

354 

337.8 

3.58 

4.68 

3.97 

DTRAN 

502 

561 

529.5 

10.71 

13.75 

11.91 

ASSIGN 

* 

* 

* 

5.02 

7.22 

6.30 

GNET 

756 

933 

812.8 

5.06 

6.17 

6.95 
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Table  4,2.  (continued) 


No 

Min 

. of  Pivots 
Max  Mean 

Min 

CPU  Times 
Max 

Mean 

Problem  AD-1 

DUAL 

54 

80 

69.0 

0.21 

0.50 

0.36 

DTRAN 

103 

124 

114.0 

0.48 

0.90 

0.70 

ASSIGN 

* 

* 

* 

0.18 

0.29 

0.25 

GNET 

164 

190 

177.6 

0.73 

0.86 

0.80 

Problem  AD -2 

DUAL 

147 

187 

169.8 

1.60 

3.04 

2.47 

DTRAN 

234 

261 

253.6 

3.02 

5.17 

4.56 

ASSIGN 

* 

* 

* 

1.28 

1.82 

1.49 

GNET 

392 

547 

448.6 

3.85 

5.80 

4.51 

Problem  AD- 3 

DUAL 

226 

255 

243.0 

4.87 

7.12 

5.60 

DTRAN 

393 

436 

418.5 

15.18 

21.69 

17.27 

ASSIGN 

* 

* 

* 

2.64 

3.82 

3.25 

GNET 

628 

714 

671.3 

9.60 

11.05 

10.26 

Problem  AD-4 

DUAL 

332 

351 

338.8 

10.99 

16.76 

14.44 

DTRAN 

519 

559 

540.3 

33.12 

52.46 

39.90 

ASSIGN 

* 

* 

* 

5.85 

9.83 

7.89 

GNET 

900 

1114 

993.3 

21.32 

23.10 

22.60 

When  we  examine  Table  4.2,  we  see  that  for  sparse  problems, 
algorithm  DUAL  is  faster  than  all  the  other  algorithms.  It  is  2-3  times 
faster  than  the  other  dual  simplex  algorithm  DTRAN,  and  1.5-2  times 
faster  than  GNET  and  ASSIGN. 

However,  for  dense  problems,  algorithm  ASSIGN  is  the  fastest 
algorithm.  It  is  1.5-2  times  faster  than  DUAL.  But,  DUAL  is  still  2-3 


times  faster  than  GNET  and  DTRAN. 


CHAPTER  5 

THE  PROJECT  SCHEDULING  PROBLEM 


Large  complex  projects  with  interdependent  activities  can  be 
described  by  project  networks.  In  a project  network  each  arc  represents 
an  activity  and  each  node  represents  an  event  signalling  the  completion 
of  all  activities  leading  into  it.  The  structure  of  the  network  defines 
the  relations  between  activities  and  events.  A schedule  associates  an 
occurrence  time  with  each  event.  Therefore,  the  project  can  be 
scheduled  in  several  different  ways.  We  assume  that  known  amounts  of 
cash  flows  are  associated  with  each  activity  and  with  some  events  in 
the  project.  Hence  given  any  feasible  schedule  the  present  value  of  all 
cash  flows  can  be  calculated.  In  this  chapter,  we  consider  a project 
scheduling  problem  introduced  by  Russel  [45] . Given  a project  network 
and  the  cash  flows  associated  with  the  activities  and  events,  the 
problem  is  one  of  scheduling  activities  and  events  to  maximize  the 
present  value  of  the  outlays  and  receipts  that  occur  during  the 
execution  of  the  project. 

5.1  Formulation  of  the  Problem 

In  a deterministic  activity  network  G with  arc  set  E and  node  set 
V,  each  arc  (i,j)eE  represents  an  activity  and  each  node  keV  represents 
an  event  signalling  the  completion  of  all  activities  leading  into  it. 
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Let  /3(i)  represent  the  "backward  star"  of  node  i which  is  all  the 
activities  leading  into  node  ieV,  and  let  °c(i)  represent  the  "forward 
star"  of  node  i which  is  all  the  activities  emanating  from  node  ieV. 
There  are  two  special  nodes  in  V;  node  1 and  node  n that  represent  the 
start  and  end  of  project,  respectively.  In  this  chapter,  it  is  assumed 
that  the  activities  are  indivisible,  i.e.,  once  an  activity  starts  it 
must  be  completed  without  any  interruption.  Moreover,  we  will  assume 
that  each  node  is  on  a directed  path  from  node  1 to  node  n.  Given  a 
start  time  of  activity  (i,j)eE  there  is  an  associated  cash  flow  Cjjt 

to  occur  at  time  T^+t,  t=0,l djj  where  djj  is  the  duration  of 

activity  (i,j).  If  C£jt  > 0 then  it  is  the  revenue,  and  if  c-jjt  < 0, 
then  it  is  the  cost  associated  with  the  activity.  Given  the  start  time 
of  an  activity  (i,j)  the  net  present  value,  NP^j  of  (i,j)  at  time 
is  given  by 

dij 

NPjj  - 2 cijt  cc-t  (5.1) 

t=0 

Here  oc'^  is  called  the  one-period  discount  factor  and  given  by  (1+i)'^ 
for  discrete  compounding  and  e'r  for  continuous  compounding,  where  i is 
called  the  interest  rate  and  r is  called  the  nominal  interest  rate. 

Similarly,  given  the  completion  time  Tj  of  an  activity  (i,j)  the 
net  future  value  of  its  cash  flows  at  time  Tj  is  given  by 


NFij  = E Cjjt  oc 
t=0 


( di j *t) 


(5.2) 


Note  that  for  each  activity  (i,j)eE,  and  NPjj  agree  in  sign.  In 

what  follows  we  assume  continuous  compounding  and  thus  oc*l  = e‘r. 
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Let  represent  the  realization  time  of  event  ieV.  Given  such  a 
feasible  realization,  for  any  (i,j)eE  we  have 

Ti  + dij  < Tj  . (5.3) 

For  an  activity  (i,j)eE,  if  NP^j  > 0,  then  clearly  (i,j)  will  be 
scheduled  to  start  at  time  T^,  otherwise  it  will  be  scheduled  to  end  at 
time  Tj . Therefore  regardless  of  the  form  of  cash  flows  on  the 
activities  we  can  construct  an  equivalent  cash  flow  problem  in  G where 
cash  flows  occur  only  at  Tj,  ieV,  where 

Fi  = Z NFki  + Z NPjj  + 9i  ieV.  (5.4) 

(k,  i)e/?(i)  (i,j)e«(i) 

NFki  < 0 NPij  > 0 

Here  ^ is  a receipt  associated  with  the  realization  of  event  ieV. 

Therefore,  for  given  F^  associated  with  node  i the  problem 
becomes  one  of  finding  realization  time  T^  of  each  event  i such  that 
the  present  value  of  all  Fj^,  ieV  is  maximized.  It  is  assumed  that  there 
is  an  upper  bound  (R)  on  the  duration  of  the  entire  project.  In  the 
case  no  such  bound  is  available  then  without  loss  of  generality  R can 
be  set  equal  to  an  arbitrary  large  positive  number. 

Let  T — (T^ , T2 , . . . , Tn) . The  project  scheduling  problem  may  be 
formulated  as  : 

n n 

(PI)  Maximize  f(T)  = Z fi(Ti)  = Z Ft  exp(-rTi)  (5.5) 

i-1  i-1 

Subject  to 


Tt  - Tj  < -d^ 

<i,j)eE 

(5.6) 

o 

II 

r— 1 

H 

(5.7) 

Tn  < R 

(5.8) 
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where  r is  the  nominal  interest  rate  used  for  discounting.  Without  loss 
of  generality,  we  have  assumed  that  the  project  starts  at  time  0 and 
must  be  completed  by  time  R.  Note  that  fi(T^)  is  a strictly  increasing 
(decreasing)  concave  (convex)  function  for  < 0 (F^  > 0)  for  T^  > 0. 

The  statement  of  the  project  scheduling  problem  as  such  leads  to 
a very  interesting  and  critical  conclusion  that  the  usual  "time 
critical  path"  may  be  irrelevant  to  the  maximizing  of  present  value; 
and  that  it  may  be  the  optimal  policy  never  to  finish  the  project.  To 
illustrate  the  situation  when  the  optimal  policy  is  never  to  finish  the 
project  we  refer  to  the  case  of  "unbalanced"  or  "front  loaded"  bidding 
as  described  in  [45]. 

A proposal  for  many  types  of  contracted  projects  includes  a list 
of  jobs  to  be  completed.  The  project  contractor  then  bids  a price  for 
each  job.  The  contract  is  awarded  to  the  contractor  whose  total  bid  is 
the  smallest.  The  contractor  is  then  paid  the  price  of  each  job  as  each 
job  is  completed.  Given  total  project  bid,  it  is  to  the  contractor's 
advantage  to  increase  the  price  for  jobs  to  be  completed  early  and 
decrease  the  price  of  those  to  be  completed  late  in  the  project  without 
changing  the  total  bid.  Not  only  does  this  strategy  increase  the 
present  value  of  the  project,  but  it  may  also  provide  an  incentive  for 
the  contractor  never  to  finish  the  project  or  delay  it  as  long  as 
possible . 

However,  if  the  contractee  was  to  use  the  net  present  value  of 
cash  flows  as  a criterion  such  an  occurrence  could  be  avoided  by 
imposing  a sufficiently  large  penalty  for  delaying  the  project.  We  also 
note  that  this  example  illustrates  a problem  which  is  not  uncommon 
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especially  in  highly  inflationary  economies.  Potential  applications  of 
the  model  are  almost  as  numerous  as  are  the  applications  of  critical 
path  methods.  Some  specific  applications  are  given  in  [45]  as 

(a)  the  pricing  of  delays  in  contractual  situations, 

(b)  the  optimum  timing  of  logically  dependent  investment 

(c)  the  avoidance  of  payment  structures  which  positively 
discourage  the  early  completion  of  a project, 

(d)  identifying  "cost  control"activities  during  the  progress  of  a 
contract , 

(e)  the  optimum  timing  of  major  projects  in  a developing  national 
economy . 


5,2  Literature  Review 

The  project  scheduling  problem  was  introduced  by  Russel  [45] . His 
approach  for  solving  problem  PI  is  to  linearize  f(T)  by  taking  Taylor 
series  expansion  of  f(T)  around  a feasible  solution  T and  taking  only 
the  linear  terms.  The  dual  of  the  linearized  problem  possesses  a 
minimum  cost  network  flow  structure.  The  dual  problem  is  solved 
iteratively  for  each  new  approximation  until  the  convergence  to  a local 
optimum  is  realized. 

Grinold  [46]  proposed  the  following  transformation  of  the 


problem: 
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n 


(P2)  Maximize  Z y^  F^ 

i-1 

(5.9) 

Subject  to: 

yj  kij 

< i 

, (l,j)eE 

(5.10) 

yi  Kij 

- yj  < ° 

, i+1,  (i,j)eE 

(5.11) 

-yn 

* -Knl 

(5.12) 

where  Kjj  — exp  (rd^j),  y^  = exp  (-rTj_).  He  further  proved  that  PI  and 
P2  are  equivalent.  There  is  a one-to-one  correspondence  between  the 
extreme  points  of  PI  and  P2 , and  for  given  T,  if  y(T)  is  optimal  to  P2 
then  T is  optimal  to  PI. 

Doersch  and  Patterson  [47]  solved  the  problem  of  maximizing  the 
present  value  of  cash  flows  subject  to  capital  constraints. 

5,3  The  Project  Scheduling  Algorithm 

The  algorithm  presented  here  is  a modification  of  the  algorithm 
developed  by  Erenguc  and  Tufekci  [48], 

5.3,1  The  Extreme  Point  Property  of  Project  Scheduling  Problem  f 48 1 

In  this  section  we  characterize  the  extreme  points  of  S and  then 
present  the  theorem  which  establishes  the  basis  for  the  proposed 
algorithm.  For  convenience  we  will  formulate  PI  as  follows: 

n 

Maximize  f(T)  = Z Ft  e'rTi  (5.13) 

i-1 

Subject  to 

TeS  = { T:  - Tj  < -dqj  for  (l.j)eE;  Tt  - Tj  < -djj  for 

i+l» (i,j)eE;  Tn  < R } 


(5.14-5.16) 
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Here  we  assume  d]j  > 0,  (l,j)eE,  and  F]^  - 0;  if  not  a dummy  node 
0 and  a dummy  activity  (0,1)  with  dQi=  e > 0 and  Fq  = 0 can  be  appended 
to  G that  will  insure  Tj  > 0,  j— 1, . . . ,n.  The  optimal  solution  to  the 
modified  problem  is  given  as  e"r£f(T*)  where  f(T*)  is  the  optimal 
solution  to  PI.  Thus,  without  loss  of  generality  we  assume  d^j  > 0,  for 
all  (1, j)eE.  In  what  follows  we  assume  |E|  = m.  We  note  that  T^  = 0 
does  not  appear  in  this  formulation.  Using  the  matrix  notation  we  have 
maximize  f(T)  (5.17) 

s.t.  AT  < d'  (5.18) 


where  A is  an  (m+l)x(n-l)  constraint  matrix,  T is  the  column  vector  of 

variables  T2 Tn,  and  d'  = [ -d,  R ]^.  Adding  slack  variables  to 

(5.18)  we  get 

Maximize  f(T)  (5.19) 

s.t.  AT  + IS  = d'  (5.20) 


We  note  that  the  matrix  A = [A, I]  has  full  row  rank.  Moreover,  any 

basic  feasible  solution  to  (5.20)  will  have  T2 Tn  basic  since 

T^  > 0,i=2,...,n.  Let  B be  any  feasible  basis  matrix  of  A.  Also  let  the 
rows  and  the  columns  of  B be  arranged  in  such  a way  that  the  first 

(n-1)  columns  and  rows  correspond  to  Tj , j=2 n.  Since  B is  an 

(m+l)x(m+l)  matrix,  there  will  be  m-n+2  slack  variables  Sg  in  the 
basis.  Thus  B can  be  put  in  the  form  given  below. 


B 


T 


SB 


T SB 


n- 1 


m-n+2 


p 

0 

Q 

1 

n- 1 


m-n+2 
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We  note  here  that  each  row  of  P corresponds  to  an  arc  (i,j)  of  the 
node-arc  incidence  matrix  of  the  project  network.  Since  B is  a basis 
matrix,  B-^  exists  which  in  turn  implies  P'^  exists.  From  the  theory 
of  network  flows  [49],  the  matrix  P induces  a spanning  tree.  After 
rearranging  and  partioning  the  right  hand  side,  the  basic  solution  can 
be  obtained  by  solving  the  following  system  in  two  stages. 


p 

0 

Q 

I 

/ 

d'p 

dB 


(5.21) 


First  by  solving 

PT  - dT'  (5.22) 

as  in  dual  variable  computation  of  network  simplex  [49] , and  then  using 
(SB)ij  = Tj  -Ti  - dij  (5.23) 

for  all  those  arcs  which  are  not  in  the  spanning  tree.  We  have  just 
shown  that  an  extreme  point  solution  to  PI  can  be  characterized  by  a 
spanning  tree  of  the  project  network.  Naturally,  for  an  extreme  point 
solution  to  be  feasible  we  must  have  Tj;  > 0,  i=2,...,n  and  Sjj  > 0 for 
all  (i,j)eE. 

In  the  remainder  of  this  section  we  will  show  that  an  optimal 
solution  to  the  network  scheduling  problem  occurs  at  an  extreme  point. 

Grinold  has  shown  that  the  optimal  solution  to  PI  (P2)  occurs  at 
an  extreme  point  of  the  polyhedral  set  S.  He  uses  the  equivalence  of  PI 
and  P2  to  prove  this  result.  In  this  section,  we  give  an  alternative 
proof  for  this  property  of  PI  which  unlike  Grinold' s proof  does  not 
depend  on  the  equivalence  of  PI  and  P2 . Our  proof  is  solely  based  on 
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the  structure  of  PI.  We  chose  to  present  this  proof  because  it  is  in 
the  spirit  of  and  lays  the  groundwork  for  the  algorithm  to  be  presented 
in  the  next  section. 

Definition  5.1.  Let  V Q V be  a node  set  with  the  following 
properties.  For  each  pair  of  nodes  i,jeV  there  exists  at  least  one 
chain  (arc  directions  ignored)  Pjj  between  the  nodes  i and  j such  that 
for  each  arc  (x,y)eP^j , (x,y)eE  we  have  Tx  - Ty  = -dXy,  and  for  each 
arc  (w,z)eE,  weV  and  z^V  or  w^V  and  zeV  we  have  Tw  - Ty  < -dwz.  We 
call  the  subnetwork  G induced  by  the  node  set  V and  the  arc  set 
E - {(i,j):  (i,j)eE,  Tj_  - Tj  - - j ) a critically  connected  component 
of  G. 

Proposition  5.1.  Given  a critically  connected  component 

» » » t t ! / 

G = (V  ,E  ),  there  exists  a tree  with  node  set  V and  arc  set  E C E 
which  spans  V such  that  for  each  arc  (iIj)€E,’l  we  have 

Ti  * Tj  - 'dij- 

Theorem  5,1.  There  exists  a global  optimum  solution  to  PI  which 
is  an  extreme  point  of  S. 

Proof : Let  TeS  be  an  optimum  solution  to  PI.  Let  Vj_,V2 V^, 

k < n be  the  critically  connected  components  corresponding  to  this 

k 

solution.  We  have  n Vi  - 4>  for  all  i=f=j  , and  u = V. 

i=l 

If  k=l,  then  V^=V  and  thus  there  exists  a spanning  tree  as  described  in 
Proposition  5.1.  Thus  choosing  a spanning  tree  from  V would 
characterize  an  extreme  point  solution. 


If  k > 1,  then  consider  any  component  V^.  Define 


“(Vi)  - { (x,y) : x^Vi(  yeVit  (x,y)eE) 

and 


/9(Vi)  - { (x,y) : xeVi,  y^Vi(  (x,y)eE) 

For  each  (x,y)eoc(V^)  we  have  Tx  - Ty  < -dXy  for  otherwise  node  x would 
have  been  in  V^.  Similarly,  for  each  (x,y)e^9(V^)  we  have 


Tx  - Ty  < -dxy. 

Let  = S Fx  e"rTx  Now  consider  the  following  solution  T* : 
xeVi 


T 


x 


Tx  ^ x#Vi 


Tx  + A if  xeV^ 


where  A is  defined  as 


(5.24) 


A 


min 

(x,y)G^(Vi) 


-min 

(x,y)ex(Vi) 


(SXy)  Spq 

if 

Pi  < 0 

(sxy)  — "Spq 

if 

Pi  > 0 

(5.25) 


Note  that  minimization  over  a null  set  is  assumed  to  yield  A = 0.  Here 
sXy  is  the  slack  on  arc  (x,y).  Let  the  minimum  occur  at  arc  (p,q).  The 
value  of  the  objective  function  at  t'  is 


f(T')  - 2 Fxe ' rTx  + 2 Fxe'r(Tx+A) 
x^V^  xgV^ 


(5.26) 


112 


Evaluating  f(T  ) - f(T)  we  get 

f(T')  - f (T)  - (e"rA  - 1)  Pi.  (5.28) 

From  the  definition  of  A,  we  have  A > 0 for  P^  < 0 and  thus 
(e'rA  - l)Pi  > 0.  Similarly,  we  have  A < 0 for  Pi  > 0 and  again 
(e'rA  - l)Pi  > 0.  Therefore  T*  is  at  least  as  good  a solution  as  T. 
Moreover  T has  at  least  one  less  critically  connected  component  since 
at  T*  the  component  Vw  where  node  q (or  p)  belongs  is  now  critically 
connected  to  node  set  Vi.  Continuing  this  argument  we  conclude  that 
there  exists  a solution  T which  is  at  least  as  good  as  T and  forms  a 
single  critically  connected  component.  Thus  any  spanning  tree  of  this 
component  as  described  in  Proposition  5.1  will  define  an  extreme 
point . 

Theorem  5,2.  Every  local  maximum  to  (PI)  is  a global  maximum. 

Proof:  See  Grinold  [46] . 

Theorems  5.1  and  5.2  present  an  interesting  result  in  that 
although  f(T)  is  not  convex,  the  global  maximum  occurs  at  an  extreme 
point  of  S and  moreover,  a local  optimum  is  a global  optimum. 

5,3,2  Notation 

Throughout  Section  5 . 3 we  will  use  the  following  notation. 

A path  in  G is  an  alternating  sequence  Z=vq , e^ , . . . , e^, v^  of  nodes 
and  arcs  such  that  e^  = (v^.^.v^)  or  e^  = (v^.v^.^).  In  the  former 
case,  e^  is  called  a forward  arc  of  Z.  In  the  latter  case  it  is  called 
a backward  arc  of  Z. 

Let  T-(V,E)  be  a spanning  tree  rooted  at  node  1.  Define  PR(x)  to 
be  the  predecessor  function  indicating  the  node  which  precedes  node  x 
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on  the  unique  path  from  node  xeV  to  node  0 in  T.  Let  D(i)  represent  a 
depth  function  defined  as  the  number  of  nodes  traversed  from  node  0 to 


node  i on  a unique  path  in  T. 

Let  Vp  be  the  set  of  nodes  ieV  such  that  the  arc  (PR(i),i)  is  a 
forward  arc  in  T.  Similarly,  let  Vg  be  the  set  of  nodes  ieV  such  that 
the  arc  (PR(i),i)  is  a backward  arc  in  T. 

5.3,3  The  Project  Scheduling  Algorithm.  PUSH 

In  this  section,  we  will  present  an  algorithm  which  will  solve 
(PI)  optimally.  Before  giving  the  precise  steps  of  the  algorithm,  we 
will  describe  the  logic  of  the  algorithm. 

Let  T - (V,E)  be  any  spanning  tree  of  the  project  network  rooted 
at  node  1.  We  will  call  a subtree  - (V^.E^),  the  subtree  "induced" 
by  node  ieVp  if  this  subtree  is  obtained  by  dropping  arc  (PR(i),i).  We 
define  P^  as  the  total  present  value  of  the  cash  flows  occuring  at 
nodes  jel^, 

Pi  - 2 F ^ e ’ rT J (5.29) 

jeVi 

where  Vi  is  the  node  set  of  the  subtree  induced  by  node  i.  The  logic  of 
the  algorithm  is  to  delay  all  the  events  in  the  induced  subtree  Ti  if 
the  present  value  of  cash  flows  Pi,  in  Ti,  is  negative. 

Given  an  induced  subtree  Ti  = (Vi,Ei),  we  define  a subtree 
Ti”  (Vi,Ei)  of  Ti  by  the  following  procedure,  PRUNE4. 
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PROCEDURE  " PRUNE4" 

INPUT  : rt  - (Vi.Ei),  D(k) , keVi  n VF. 

OUTPUT  : A subtree  (V^,E^)  Q with  the  property  that  for  each 

jeVip  Pj  > 0,  and  its  associated  Pi  value. 

BEGIN 

INITIALIZE  : Set  < rit  R=(j=fi6VinVF:Pj<0), 

WHILE  R + 4>  DO 

Find  jeR  such  that  D(j)  > D(k)  , for  all  k=|=jeR.  Obtain 
rr  (Vj  - Ej  > • Set  ri  - (Vi  I Vj  , E-|Ej  u { (PR(  j ) , j ) ) )■ 
Update  Pj , jeVFn  VF  and  R. 

ENDWHILE 

Compute  p£  - E Fje’rTj  (5.30) 

jevi 

END 

Here  P^  might  be  viewed  as  the  individual  contribution  of  subtree 
to  the  negative  present  value  PF.  The  algorithm  advances  the  nodes  in 
FF  if  PF  <0,  that  is  the  individual  contribution  of  rj.  is  negative. 

If  Pi  < 0 but  Pi  >0,  then  the  negativeness  is  caused  by  a set  of 
nodes  in  an  induced  subtree  rk  where  ^ C Ti,  and  hence  only  the  nodes 
in  rk  should  be  advanced. 

Let  Q be  the  set  of  nodes  with  Pi  <0.  The  algorithm  starts  with 
the  longest  path  tree  (early  start  tree)  r®,  and  finds  a node  peQ  with 
Pp  <0  such  that  no  other  node  j on  the  unique  path  from  node  p to 
node  1 will  have  Pj  <0.  Initially,  an  easy  way  of  determining  Tp  , peQ 
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is  to  check  the  depth  functions  of  the  nodes  in  Q and  select  Tp  such 

that  D(p)-  min  D(j). 

jeQ 

Throughout  the  algorithm,  after  a particular  pivot,  k,  Tp  is 
determined  while  updating  . Let  Wp»  be  the  set  of  nodes  on  the 

unique  path  from  node  p found  in  pivot  k to  node  1 in  r^+^.  Then  new  p 
is  chosen  to  be  the  node  oceWp>nQ  with  the  smallest  depth  function.  If 
Wp-nQ  = <j> , then  p is  chosen  to  be  the  node  oceQ  with  the  smallest  depth 
function.  If  no  such  node  p exists,  then  the  algorithm  terminates. 
Otherwise,  the  Tx  values  of  the  nodes  xeVp  are  advanced,  i.e.,  the 
events  are  delayed,  by  A > 0. 

The  quantity  A is  determined  by  the  slacks  of  the  out  of  tree 
arcs  from  node  set  Vp  to  node  set  V\Vp . If  no  such  arc  exists,  then  A 
is  set  to  A=R-Tn  where  R is  the  latest  delivery  time  of  the  project. 

If  there  exists  some  arcs  from  node  set  Vp  to  node  set  V\Vp,  A is 
obtained  from 

A = min  (Ty  - Tx  - dxy  : xeVp , yeV|Vp,  (x,y)eE) . (5.31) 

jhe  algorithm  obtains  the  adjacent  extreme  point  by  dropping  the  arc 
(PR(p) , p)  out  of  T and  adding  the  arc  (n,l)  if  A=R-Tn  or  adding  the  arc 
(r,s)  which  defines  A in  equation  (5.31),  into  the  tree,  T.  We  now  give 
a formal  statement  of  the  algorithm. 


ALGORITHM  PUSH 

INPUT:  G - (V,E),  F^,  ieV,  djj  , (i,j)eE  and  R,  the  project  deadline, 

r,  the  nominal  interest  rate. 

OUTPUT:  Optimal  timings  T^*  of  events  ieV  and  optimal  cost  f(T*). 
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PROCEDURE  : 
INITIALIZE: 


MAIN  STEP  : 


Find  a longest  path  tree  T®  of  G and  let  T®  be  the  vector 
of  path  lengths  from  node  1 to  all  other  nodes. 

I It 

For  each  node  jeVp  determine  the  subtree  T j - (Vj , Ej ) 
compute  Pj  by  Procedure  PRUNE4 . 

Define  Q - { j : P j < 0,  jGVp} . Set  p’=0. 

WHILE  Q=f  4>  DO 
BEGIN 

IF  p'-O  THEN 

Find  a node  peQ  such  that  D(p)  = min  { D ( j ) } . 

j^Q 

ELSE 

Find  a node  peQnWp»  such  that  D(p)=  min  { D ( j ) } 

jeQnWp, 

ENDIF 

Define  X - {(i,j):  ieVp,  j£vp>  (i,j)eE) 

IF  X - <f>,  THEN 
A = R - Tn. 

Set  T^  = Ti  + A ieVp 

Ti  - Tp  ieV\Vp . 

Update  T by  deleting  arc  (PR(p),p)  and  adding  arc 
(n, 1) . 

Update  Pj , Q. 

ELSE 

Compute 

A - min  ( Tj  - T^  - d^j  : (i,j)eX).  Let  the  minimum 
occur  at  (r,s)eX. 
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Set  Tj_  - + A ieVp 

Ti  - Ti  ieV\Vp . 

Update  T by  deleting  arc  (PR(p),p)  and  adding  arc 
(r , s) . 

Update  Pj , Q. 

END  IF 

IF  WpnQ-  <t>,  THEN 
Set  p'-O 

ELSE 

Set  p'=p 

END  IF 
ENDWHILE 

Compute  f(T)  and  print  T and  f(T). 

END. 


We  call  each  execution  of  the  main  step  of  this  algorithm  a 
"PUSH",  since  we  are  pushing  the  nodes  of  the  subtree  Tp  always  in  one 
direction.  We  note  that  each  execution  of  the  main  step  corresponds  to 
a simplex  pivot  in  the  sense  that  the  algorithm  moves  to  an  adjacent 
extreme  point  in  S as  a result  of  the  "PUSH"  operation. 

Example  5.1.  The  project  scheduling  problem  with  6 nodes  is 
solved  in  Figure  5.1.  The  nominal  interest  rate,  r,  is  given  as  .10, 
and  the  project  deadline,  R,  is  given  as  20. 
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-20 


T=(T1)T2>T3,T4,T5lT6)  = (0,4,6,3,10,12) 

P2  - -20e--1(4)  - -13.4064 

P3  . -60e-  .1(6)  - 30e'  • 1^10)  + ZOe*-1^2)  - -37.9412 
P4  - 80  e" • - 59.2654 
P5  = -30e*  • -L^10)  + 20e~  • 1(12)  - -5.0125 
P6  = 20e' • 1(12)  - 6.0239. 


Figure  5.1  An  Example 
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p6  ” p6-  p5=p5  > p4=p4  ■ p2  = p2- 
P3  - -60e- .1(6)  - -32.9287. 

Q - { 2,3,5). 

Iteration  1.  p - 2 , arc  (1,2)  leaves,  arc  (2,6)  enters,  A = 1. 


T = (0,5,6,3,10,12) 

P6  = -20e* • 1(5)  + 20e' • p(12)  = -6.1067 

P5  - -20e- . 1(5)  + 20e*-1(12)  - lOe'-1^10)  - -17.1431 

P3  - -20e- .1(5)  + 20e--1(12)  - SOe’-1^10)  - eOe'-1^6)  = -50.0687. 

P4  = 59.2654. 

p6  = p6-  p4  = p4> 

P5  = -11.0333 
P3  = -32.9287. 

Iteration  2.  p=3 , arc  (1,3)  leaves.  X=</>,  A = R-Tn  = 20-12=8,  arc  (6,1) 
enters . 


Figure  5.1  (continued) 
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T = (0,13,14,3,18,20) . 

=-  59.2654.  Q = <f> . Current  solution  T is  optimal. 

Figure  5.1.  (continued) 


5.4,  Computational  Experience  for  Project  Scheduling  Problem 

In  this  section,  we  present  the  computational  experience  obtained 
by  the  algorithm  "PUSH"  and  compare  it  with  Grinold's  algorithm. 

The  difference  between  the  algorithm  "PUSH"  and  Grinold's 
algorithm  is  the  selection  of  the  leaving  arc.  In  Grinold's  algorithm 
an  arc  (x,y)  is  selected  as  a leaving  arc  if  wXy  < 0.  Here  wjj  is 
computed  by  solving  the  following  set  of  equations: 

S w j £ - E w^j  = b^  for  i=2,3 n (5.32) 

je/3i  j eri 

where  b^  = Fj_exp( ‘rTi)  , is  the  set  of  events  immediately  preceding 
event  i and  is  the  events  immediately  following  event  i. 

Since  no  selection  rule  is  indicated  in  Grinold's  paper,  we 
assumed  two  different  selection  rules,  one  being  the  most  negative  wjj  , 
and  the  other  being  the  first  found  negative  wjj  . We  coded  two 
algorithms  using  these  selection  rules  and  named  them  GRINOLD1,  which 
implements  the  first  selection  rule  (most  negative  Wjj ) , and  GRIN0LD2 , 
which  implements  the  second  selection  rule  (first  encountered  negative 
wij). 

Algorithm  "PUSH"  chooses  an  arc  (x,y)  to  leave  the  basis  if  (x,y) 
defines  an  induced  subtree  Tp  with  Pp  < 0,  and  D(p)=  min  (D(j),  jeQ)  . 
Here  Tp  does  not  contain  any  induced  subtree  Tj  C Tp  with  Pj  <0.  This 
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is  equivalent  to  finding  an  arc  (x,y)  with  a negative  wXy  found  by 
equations  (5.32)  after  pruning  negative  j ' s computed  earlier. 

All  of  the  algorithms  tested  are  coded  in  FORTRAN  77  and  are  run 
on  a VAX  11/750  using  the  VAX-11  Fortran  V3.0  compiler.  All  runs  are 
made  overnight  when  we  were  the  only  users  of  the  system. 

We  have  tested  the  three  algorithms  for  12  sets  of  test  problems 
for  n e {50,100,150,200}  and  densities  9 - {0.1, 0.3, 0.5),  and  used  10 
replications  in  each  set.  Here  9 - 2|E|/n(n-l)  is  the  ratio  of  the 
number  of  actual  arcs  in  the  network  to  maximum  possible  number  of 
arcs.  In  all  problems  a nominal  interest  rate  of  r-0 . 10  is  used  , and 
cash  flows,  F-l,  are  generated  uniformly  between  (-10000,10000)  for  i=}=  n 
and  Fn  is  set  to  30,000. 

The  characteristics  of  the  problems  are  summarized  in  Table  5.1. 
The  performance  of  the  algorithms  are  presented  in  Table  5.2  where  the 
CPU  times  are  in  seconds,  and  the  mean  values  are  obtained  from 
Table  5.1.  Characteristics  of  the  Problems  Generated 


Problem  Type 

Number  of  Nodes 

Number  of  Arcs 

P-1 

50 

123 

P-2 

50 

369 

P-3 

50 

615 

P-4 

100 

495 

P-5 

100 

1485 

P-6 

100 

2475 

P-7 

150 

1118 

P-8 

150 

3353 

P-9 

150 

5590 

P-10 

200 

1990 

P-11 

200 

5970 

P-12 

200 

9950 
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Table  5.2.  Performance  of  the  Algorithms 


1 

No. of  Pivots 

CPU  times 

Min 

Max 

Mean 

Min 

Max 

Mean 

Problem  P-1 

PUSH 

21 

29 

25.4 

0.20 

0.38 

0.31 

GRIN0LD1 

22 

36 

26.2 

1.10 

1.77 

1.29 

GRIN0LD2 

44 

100 

60.8 

1.79 

3.91 

2.51 

Problem  P-2 

PUSH 

35 

44 

41.3 

0.51 

0.81 

0.63 

GRIN0LD1 

38 

53 

44.8 

1.98 

2.62 

2.26 

GRIN0LD2 

62 

255 

170.0 

2.65 

11.14 

7.20 

Problem  P-3 

PUSH 

35 

49 

43.8 

0.62 

0.98 

0.80 

GRIN0LD1 

39 

54 

47.1 

1.76 

3.05 

2.46 

GRIN0LD2 

84 

312 

194.3 

3.74 

13.52 

8.51 

Problem  P-4 

PUSH 

55 

99 

78.2 

1.50 

3.29 

2.44 

GRIN0LD1 

61 

104 

83.5 

10.23 

17.42 

13.98 

GRIN0LD2 

181 

332 

234.9 

23.52 

49.36 

35.48 

Problem  P-5 

PUSH 

82 

109 

97.4 

2.90 

4.73 

4.00 

GRIN0LD1 

94 

123 

109.0 

16.68 

21.22 

19.15 

GRINOLD2 

399 

782 

678.0 

64.47 

123.89 

103.54 

Problem  P-6 

PUSH 

74 

116 

93.6 

3.51 

6.82 

4.81 

GRIN0LD1 

85 

123 

110.3 

16.81 

22.53 

20.12 

GRINOLD2 

499 

1484 

936.8 

78.20 

218.79 

120.82 

Problem  P-7 

PUSH 

132 

161 

142.4 

7.00 

12.83 

9.43 

GRINOLD1 

137 

193 

163.6 

60.88 

77.96 

70.29 

GRINOLD2 

416 

887 

575.7 

161.58 

324.56 

220.04 

123 


Table  5.2.  (continued) 


Min 

No. of  Pivots 

Max  Mean 

Min 

CPU  times 
Max 

Mean 

Problem  P-8 

PUSH 

125 

174 

151.2 

9.41 

14.92 

12.28 

GRIN0LD1 

167 

228 

187.3 

69.49 

94.86 

78.00 

GRIN0LD2 

384 

1988 

1371.6 

340.20 

420.86 

509.10 

Problem  P-9 

PUSH 

137 

206 

161.4 

11.74 

20.34 

15.42 

GRIN0LD1 

158 

268 

196.8 

70.73 

117.65 

90.83 

GRIN0LD2 

1729 

4026 

2701.0 

664.31 

1482.73 

944.13 

Problem  P-10 

PUSH 

129 

229 

175.6 

9.24 

25.06 

17.23 

GRINOLD1 

205 

285 

226.3 

147.60 

201.25 

175.37 

GRIN0LD2 

488 

1106 

858.5 

364.54 

897.70 

645.97 

Problem  P-11 

PUSH 

168 

243 

197.7 

21.23 

34.98 

23.39 

GRIN0LD1 

216 

297 

266.1 

160.38 

218.80 

205.28 

GRIN0LD2 

1859 

2527 

2148.8 

1204.86 

1652.33 

1395.35 

Problem  P-12 

PUSH 

166 

211 

186.8 

24.94 

36.79 

31.16 

GRIN0LD1 

193 

348 

257.0 

148.47 

263.27 

195.78 

GRIN0LD2 

1025 

5300 

3880.1 

705.94 

4341.96 

2374.40 

When  we  examine  Table  5.2,  we  see  that  algorithm  PUSH  is  much 
faster  than  the  other  two  algorithms.  It  is  3-10  times  faster  than 
GRIN0LD1  and  8-76  times  faster  than  GRIN0LD2 . 

One  of  the  reasons  for  having  high  CPU  times  for  algorithms 
GRIN0LD1  and  GRIN0LD2  is  that  in  these  algorithms  all  w^j  ' s are 

/ 

recalculated  after  each  pivot  whereas  in  algorithm  PUSH  only  those  Pj 
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values  that  are  affected  by  the  pivot  are  recalculated.  Moreover,  in 
algorithm  GRIN0LD1,  for  determining  the  leaving  arc,  all  the  w-jj  values 
are  checked  to  find  the  most  negative  one  and  this  increases  the 
computational  time. 

The  reason  for  having  larger  number  of  pivots  in  algorithm  GRIN0LD1 
is  due  to  the  selection  rule  of  min  {w^j } . An  arc  (i,j)  may  have  the  most 
negative  wjj  value,  however,  the  negativeness  can  be  caused  because  of 
other  arcs.  In  our  terms,  Pj  may  be  positive  but  Pj  < 0 due  to  the 
P^  values  of  subtrees  I\  C Tj  with  P^  < 0.  In  such  a case,  the 
algorithm  requires  a "pulling  back"  pivot  where  events  in  Tj  are  advanced 
(not  delayed).  Algorithm  GRIN0LD2  requires  no  "pulling  back"  pivots, 
however  the  number  of  pivots  required  is  enourmously  large  compared  to 
the  other  two  algorithms.  The  reason  for  this  is  that  the  algorithm  can 
"push"  separately  several  subtrees  that  are  on  a given  path  if  each  of 
these  subtrees  have  negative  P^  values,  while  algorithm  PUSH  can  "push" 
all  of  these  subtrees  with  a single  pivot  . 


CHAPTER  6 
SUMMARY 


In  this  dissertation  we  focused  on  some  special  cases  of  optimal 
differential  problems. 

In  Chapters  2-5,  we  developed  three  dual  simplex  algorithms  for 
the  minimum  cost  flow  problem,  and  its  special  cases,  the 
transportation  problem  and  the  assignment  problem. 

The  leaving  arc  selection  rule  in  the  dual  simplex  algorithm, 
DUAL,  provided  us  to  find  a computational  bound  for  the  algorithm  of 
0 (M | V | 2 ) where  M is  the  sum  of  the  supplies  (or  demands)  in  the 
problem.  By  applying  the  scaling  procedure  developed  by  Edmonds  and 
Karp  [28],  we  developed  the  dual  simplex  algorithm, DUALSC  which 
provided  a polynomial  computational  bound  of  0[(r+l)|Vp].  Since  this 
bound  is  dependent  on  the  logarithm  of  the  input  data  of  supply  and 
demand  values,  DUALSC  is  not  a genuinely  polynomial  algorithm.  However, 
in  terms  of  the  magnitude  of  the  bound,  DUALSC  provides  the  smallest 
value  when  we  compare  it  with  the  bounds  of  other  genuinely  polynomial 
bounds  of  algorithms  developed  by  Orlin  [29]  and  polynomial  bound  of 
the  algorithm  DTRANSC,  developed  by  Ikura  and  Nemhauser  [32].  Orlin' s 
two  algorithms  have  computational  bounds  of  0( | V | ^log | V | ) and 
0( | | log | V | ) . DTRANSC  has  a computational  bound  of  0 [ (t+1) | V | ^ ] . 
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In  our  computational  experimentation,  we  compared  the  algorithms 
DUAL  and  DUALSC  with  other  dual  simplex  algorithms  and  the  primal 
simplex  algorithm,  GNET.  In  previous  experimentations  reported  in  the 
literature,  the  primal  simplex  algorithms  were  found  to  be 
computationally  superior  to  all  other  algorithms  for  large  sized 
problems.  In  our  experimentation,  for  minimum  cost  flow  problems  and 
transportation  problems,  algorithms  DUAL  and  DUALSC  performed  better 
than  all  the  other  dual  simplex  algorithms  for  both  small  and  large 
sized  problems.  They  were  3-5  times  faster  than  the  algorithms 
developed  by  Orlin,  and  1.5-2  times  faster  than  DTRANSC.  Even  though 
algorithm  DUAL  was  competitive  with  GNET  for  large  dense  problems , GNET 
performed  better  in  large  sparse  problems  by  being  2-3  times  faster 
than  DUAL. 

When  algorithm  DUAL  is  used  to  solve  assignment  problems,  it 
becomes  equivalent  to  Balinski  s [43]  dual  simplex  algorithm  and 
provides  the  same  bound  of  (N-l)(N-2)/2  pivots  where  N=|Vs|=|Vp|. 
However  Balinski  s algorithm  cannot  be  extended  to  solve  the  minimum 
cost  flow  problem  or  the  transportation  problem. 

In  our  computational  experimentation  for  assignment  problems,  we 
concluded  that  algorithm  DUAL  is  2-3  times  faster  than  other  dual 
simplex  algorithm,  DTRAN  and  the  primal  algorithm  GNET.  When  we 
compared  DUAL  s performance  with  a primal -dual  algorithm  developed  by 
Burkard  and  Derigs  [44],  we  concluded  that  DUAL  is  1.5-2  times  faster 
than  the  primal-dual  algorithm  for  sparse  problems  but  the  primal-dual 
algorithm  is  1.5-2  times  faster  than  DUAL  for  dense  problems. 
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There  still  remains  an  interesting  open  question  as  to  whether 
there  exists  a simplex  pivot  rule  that  solves  the  minimum  cost  flow 
problem  in  polynomial  time. 

In  Chapter  5,  we  developed  an  algorithm,  PUSH,  that  solves  the 
project  scheduling  problem  introduced  by  Russel  [45] . This  project 
scheduling  problem  has  an  objective  function  of  maximizing  the  present 
value  of  cash  flows . 

The  difference  between  the  algorithm  PUSH  and  the  algorithm 
developed  by  Grinold  [46]  to  solve  the  same  problem  is  the  selection 
rule  of  the  leaving  arc  in  a given  pivot.  Since  no  selection  rule  is 
indicated  in  Grinold  s algorithm,  we  assumed  two  different  commonly 
used  selection  rules  and  compared  these  algorithms  with  algorithm  PUSH 
in  our  computational  experimentation.  Based  on  our  experimentation,  we 
concluded  that  algorithm  PUSH  is  3-10  times  faster  than  the  algorithm 
which  employs  the  leaving  arc  selection  rule  of  the  most  negative 
value,  and  8-76  times  faster  than  the  algorithm  which  employs  the 
leaving  arc  selection  rule  of  the  first  found  negative  value. 


APPENDIX  A 
LIST  OF  NOTATIONS 


Chapters  2-4 


Pi 

: the  unique  path  connecting  node  i to  node  0 in  T. 

T = (V,ET) 

: a spanning  tree  rooted  at  node  0. 

Ti  - (Vi.Ei) 

: a subtree  for  node  ieVp  obtained  by  dropping 
arc  (PR(i),i)  out  of  T and  taking  the  subtree 
that  does  not  contain  the  root  node. 

Ti  - (V-,Ei)  C Ti 

: a subtree  for  node  ieVp  when  all  the  subtrees 

l f 

Tj  £ Ti  with  positive  Sj  values  are  pruned 
from  Ti. 

ft  ft  ft 

Ti  -(Vi  ,Ei  ) 

: a subtree  for  node  i€Vg  obtained  by  dropping 
arc  (i,PR(i))  out  of  T and  taking  the  subtree 
that  does  not  contain  the  root  node. 

t t t ft!  Iff  ft 

Ti  = (vi  ,Ei  ) C Ti 

: a subtree  for  node  ieVg  when  all  the  subtrees 
Tj  with  Sj  > 0 are  pruned  from  T^/  . 

Ti 

: a subtree  for  node  ieVp  when  all  the  subtrees 

f 1 • . Ilf  f 

Tj  with  nonpositive  Sj  are  pruned  from  Tp. 

Si 

: negative  of  the  flow  on  the  forward  arc  (PR(i),i) 
for  ieVp. 

1 

Si 

: negative  of  the  flow  on  the  forward  arc  (PR(i),i) 
in  Ti  for  ieVp. 

1 1 

Si  : 

: flow  on  the  backward  arc  (i,PR(i))  for  i eVg. 
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Zi 


t 

Zi 

Zi 


vF 


Vfi 


Q 

KCQ 

(.)b 


(.)a 


(.)z 


: flow  on  the  backward  arc  (i,PR(i))  for  ieVg  in 

I t I 

Ti  • 

: negative  of  the  flow  on  the  forward  arc  (PR(i),i) 
in  for  ieVp. 

: set  of  nodes  j , with  Sj  > 0 that  are  pruned  from 
subtree  or  T-j/  . 

: Zi  U i,  if  S^  > 0 ; if  < 0. 

: set  of  nodes  j with  Sj  <0  that  are  pruned  from 
subtree  T^. 

: set  of  nodes  such  that  the  arc  (PR(i),i)  is  a 
forward  arc . 

: set  of  nodes  such  that  the  arc  (PR(i),i)  is  a 
backward  arc . 

: set  of  nodes  j with  Sj  > 0. 

: set  of  nodes  j with  Sj  > 0 and  there  exists  no 
other  kePj  with  S^  > 0. 

: a value (.)  or  a subtree  (.)  before  a particular 
pivot . 

: a value(.)  or  a subtree  (.)  after  a particular 
pivot . 

: a value (.)  or  a subtree  (.)  in  the  subproblem  z. 
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